diff --git a/dp1/euclid_q1_lenses.ipynb b/dp1/euclid_q1_lenses.ipynb new file mode 100644 index 0000000..f120ebe --- /dev/null +++ b/dp1/euclid_q1_lenses.ipynb @@ -0,0 +1,926 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "f11fc55e-df06-47dd-8599-eb78bc578862", + "metadata": {}, + "source": [ + "# Euclid Q1 Lenses\n", + "\n", + "* **Phil Marshall, Phil Holloway, Ralf Kaehler, Ferro Shao**\n", + "* DP1\n", + "* data.lsst.cloud\n", + "* r29.1.1\n", + "* Fri July 11 2025" + ] + }, + { + "cell_type": "markdown", + "id": "65367d2a-7335-4752-9565-3861387a2a63", + "metadata": {}, + "source": [ + "## Goals\n", + "\n", + "* Extract _ugrizy_ coadd image cutouts for each Euclid Q1 strong lens candidate in the ECDFS and EDFS DP1 fields\n", + "* Visualize them as _gri_ color composites.\n", + "* Stretch: deconvolve them using the Rubin SharPy by Kaehler et al (in prep)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "729bff4a-1e33-458b-83c5-9921ca77fed8", + "metadata": {}, + "outputs": [], + "source": [ + "from lsst.daf.butler import Butler\n", + "import lsst.afw.display as afw_display\n", + "import lsst.geom as geom\n", + "from lsst.afw.image import MultibandExposure\n", + "import numpy as np\n", + "from astropy.visualization import make_lupton_rgb\n", + "from astropy.table import Table, Column\n", + "import matplotlib.pyplot as plt\n", + "\n", + "afw_display.setDefaultBackend('matplotlib')" + ] + }, + { + "cell_type": "markdown", + "id": "f26c3831-a5d6-4358-ba5a-4b54d06a4ed7", + "metadata": {}, + "source": [ + "## Cutout Image Extraction\n", + "\n", + "First we need to make a list (or better, a table) of targets. Then, for each one, we find out which DP1 coadd patch it lies in. (We'll need to choose which patch, for systems that lie in the patch overlap regions and hence in multiple patches.) Then, we loop over patches and bands, uploading a patch image and extracting all the cutouts we can - which will mean getting the image coordinates for each system" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ba8ee8ee-7cb4-4a9d-a3e8-186c0c593095", + "metadata": {}, + "outputs": [], + "source": [ + "butler = Butler(\"dp1\", collections=\"LSSTComCam/DP1\")\n", + "assert butler is not None" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "39bed3db-4271-4aeb-82c4-a23bab785d8e", + "metadata": {}, + "outputs": [], + "source": [ + "butler.get_dataset_type('deep_coadd').dimensions.required" + ] + }, + { + "cell_type": "markdown", + "id": "6fb6113c-2763-4424-8856-06d79b1c5de3", + "metadata": {}, + "source": [ + "### Single Sky Position, Single Band\n", + "\n", + "Let's just try extracting a single 32x32 pixel cutout image in one band." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b20bc42a-e4af-4768-ba96-6f10879314d0", + "metadata": {}, + "outputs": [], + "source": [ + "ra = 59.626134\n", + "dec = -49.06175\n", + "\n", + "band = 'i'" + ] + }, + { + "cell_type": "markdown", + "id": "248142d6-1e6c-473d-99a7-819b89d77034", + "metadata": {}, + "source": [ + "Turn the coordinates into an IAU standard object name, we'll need this later:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7dc87bd9-271c-4205-a2af-c71381cee28e", + "metadata": {}, + "outputs": [], + "source": [ + "import astropy.units as u\n", + "from astropy.coordinates import SkyCoord\n", + "\n", + "# Example RA and Dec coordinates\n", + "rahms = ra * u.hourangle # 12 hours, 30 minutes, 36 seconds\n", + "decdms = dec * u.deg # +12 degrees, 24 minutes, 0 seconds\n", + "\n", + "# Create a SkyCoord object\n", + "coordinates = SkyCoord(ra=rahms, dec=decdms, frame='icrs')\n", + "\n", + "# Format the coordinates into an IAU-style string\n", + "name = (f'EUCLID J{coordinates.ra.to_string(unit=u.hourangle, sep=\"\", precision=1, pad=True)}'\n", + " f'{coordinates.dec.to_string(sep=\"\", precision=0, alwayssign=True, pad=True)}') #\n", + "\n", + "print(name)" + ] + }, + { + "cell_type": "markdown", + "id": "7dd657db-3028-4897-ae44-a4e44c4aff63", + "metadata": {}, + "source": [ + "We need to find the tract and patch that this target is in. This approach was adopted from the CST Tutorial [\"03a Image Display and Manipulation\"](https://github.com/lsst/tutorial-notebooks/blob/main/DP0.2/03a_Image_Display_and_Manipulation.ipynb)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "34975668-16b9-447e-b37a-56a9e78cc45e", + "metadata": {}, + "outputs": [], + "source": [ + "radec = geom.SpherePoint(ra, dec, geom.degrees)\n", + "cutoutSize = geom.ExtentI(32, 32)\n", + "\n", + "skymap = butler.get(\"skyMap\")\n", + "tractInfo = skymap.findTract(radec)\n", + "patchInfo = tractInfo.findPatch(radec)\n", + "\n", + "patch = tractInfo.getSequentialPatchIndex(patchInfo)\n", + "tract = tractInfo.getId()\n", + "\n", + "dataId = {'tract': tract, 'patch': patch, 'band': band}" + ] + }, + { + "cell_type": "markdown", + "id": "1bd4189a-0cc7-4fa3-b219-ccfb57a614f2", + "metadata": {}, + "source": [ + "When testing, it can be useful to upload the whole patch image and inspect it." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "774cb05c-c831-4c79-b891-aa00db25f7be", + "metadata": {}, + "outputs": [], + "source": [ + "# deep_coadd = butler.get('deep_coadd', band=band, tract=tract, patch=patch)\n", + "# coadd" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c8f5d83d-beb8-4a01-a6c5-9a9d16a867d5", + "metadata": {}, + "outputs": [], + "source": [ + "# fig = plt.figure(figsize=(6,6))\n", + "# display = afw_display.Display(frame=fig)\n", + "# display.scale('asinh', 'zscale')\n", + "# display.mtv(coadd.image)\n", + "# plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "1d516446-ae23-46e1-bc00-0c8cdb3c6ce5", + "metadata": {}, + "source": [ + "Now to define a small bounding box, and extract the pixels in it. This first cell below _should_ work, but doesn't - maybe some tract/patch confusion. There could be some speed up here at some point, making multiple cutouts from the same patch image using the calexp object's native factory method. See the Cutout Factory demo notebook by Melissa Graham at https://github.com/lsst/cst-dev/blob/main/MLG_sandbox/random/cutout_factory_demo_2025-06-05.ipynb for a working example of making multiple cutouts from the same patch image." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e32d9519-f4ec-4499-b6bb-f15eb3bad538", + "metadata": {}, + "outputs": [], + "source": [ + "# xy = geom.PointI(tractInfo.getWcs().skyToPixel(radec))\n", + "# bbox = geom.BoxI(xy - cutoutSize // 2, cutoutSize)\n", + "\n", + "# cutout = coadd.Factory(coadd, bbox)" + ] + }, + { + "cell_type": "markdown", + "id": "bb06ca5f-af1b-4c61-8196-0690e743028f", + "metadata": {}, + "source": [ + "Here's some code that does work: define the bounding box, then just grab that part of the image." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8c09ccc4-d4a0-4091-953e-b0bec9f3e962", + "metadata": {}, + "outputs": [], + "source": [ + "xy = geom.PointI(tractInfo.getWcs().skyToPixel(radec))\n", + "bbox = geom.BoxI(xy - cutoutSize // 2, cutoutSize)\n", + "\n", + "parameters = {'bbox': bbox}\n", + "\n", + "cutout = butler.get(\"deep_coadd\", parameters=parameters, dataId=dataId)" + ] + }, + { + "cell_type": "markdown", + "id": "9b6105f8-20ad-4e6a-9000-ad468abc7978", + "metadata": {}, + "source": [ + "Quick look to check we got our object!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1dff2724-6f93-4e35-83da-f88a98d55c9e", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(3, 3))\n", + "display = afw_display.Display(frame=fig)\n", + "display.scale('asinh', 'zscale')\n", + "display.mtv(cutout.image)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "a2f16999-1a5a-45fd-a9ef-cf1714aff8fb", + "metadata": {}, + "source": [ + "### Single Object, Multiple Bands\n", + "\n", + "Loop over all 6 bands and extract the cutout image in each one." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c8c2c390-2374-41ba-a19b-c3e99c8dcd53", + "metadata": {}, + "outputs": [], + "source": [ + "bands = [\"u\",\"g\",\"r\",\"i\",\"z\",\"y\"]\n", + "cutout = {}\n", + "\n", + "for band in bands:\n", + " dataId = {'tract': tract, 'patch': patch, 'band': band}\n", + " cutout[band] = butler.get(\"deep_coadd\", parameters=parameters, dataId=dataId)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "41514463-c5c6-405f-a25c-bec43041139f", + "metadata": {}, + "outputs": [], + "source": [ + "cutout" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "533d3a49-967c-4a61-af6b-9d96a49c0bb9", + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(3, 3))\n", + "display = afw_display.Display(frame=fig)\n", + "display.scale('asinh', 'zscale')\n", + "display.mtv(cutout[\"g\"].image)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "c9b12a89-ef2c-4a62-af7f-6b8630eb899b", + "metadata": {}, + "source": [ + "OK - we have 6 cutouts for this target, so can go ahead and make a color composite - see below for this. It took about 5 secs to make them all: we'll need to keep an eye on this, and return to the `factory` approach to try and speed things up a bit." + ] + }, + { + "cell_type": "markdown", + "id": "7d216083-ef68-4d13-861d-0466311b55e5", + "metadata": {}, + "source": [ + "## _gri_ Composite Image Visualization\n", + "\n", + "Here's a useful function, adapted from the CST Tutorial [\"03a Image Display and Manipulation\"](https://github.com/lsst/tutorial-notebooks/blob/main/DP0.2/03a_Image_Display_and_Manipulation.ipynb)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6bcfa278-84c5-4eda-852d-a31e9e1f2609", + "metadata": {}, + "outputs": [], + "source": [ + "def showRGB(image, bgr=\"gri\", ax=None, fp=None, figsize=(8,8), stretch=100, Q=1, name=None):\n", + " \"\"\"Display an RGB color composite image with matplotlib.\n", + " \n", + " Parameters\n", + " ----------\n", + " image : `MultibandImage`\n", + " `MultibandImage` to display.\n", + " bgr : sequence\n", + " A 3-element sequence of filter names (i.e. keys of the exps dict) indicating what band\n", + " to use for each channel. If `image` only has three filters then this parameter is ignored\n", + " and the filters in the image are used.\n", + " ax : `matplotlib.axes.Axes`\n", + " Axis in a `matplotlib.Figure` to display the image.\n", + " If `axis` is `None` then a new figure is created.\n", + " figsize: tuple\n", + " Size of the `matplotlib.Figure` created.\n", + " If `ax` is not `None` then this parameter is ignored.\n", + " stretch: int\n", + " The linear stretch of the image.\n", + " Q: int\n", + " The Asinh softening parameter.\n", + " name: str\n", + " The name of the object/field to be displayed.\n", + " \"\"\"\n", + " # If the image only has 3 bands, reverse the order of the bands to produce the RGB image\n", + " if len(image) == 3:\n", + " bgr = image.filters\n", + " # Extract the primary image component of each Exposure with the .image property, and use .array to get a NumPy array view.\n", + " rgb = make_lupton_rgb(image_r=image[bgr[2]].array, # numpy array for the r channel\n", + " image_g=image[bgr[1]].array, # numpy array for the g channel\n", + " image_b=image[bgr[0]].array, # numpy array for the b channel\n", + " stretch=stretch, Q=Q) # parameters used to stretch and scale the pixel values\n", + " if ax is None:\n", + " fig = plt.figure(figsize=figsize)\n", + " ax = fig.add_subplot(1,1,1)\n", + " \n", + " plt.axis(\"off\")\n", + " ax.imshow(rgb, interpolation='nearest', origin='lower')\n", + "\n", + " if name is not None:\n", + " plt.text(0, 31, name, color='white', fontsize=12, horizontalalignment='left', verticalalignment='top')\n", + " \n", + " plt.text(0, 2, 'astropy Lupton '+bgr, color='white', fontsize=12, horizontalalignment='left', verticalalignment='top')\n", + " \n", + " plt.tight_layout();\n", + "\n", + " return" + ] + }, + { + "cell_type": "markdown", + "id": "a51c625e-f017-423c-9f65-753e1f7a8fb5", + "metadata": {}, + "source": [ + "First we need to package our cutouts into a MultibandExposure object, then we pass that to the RGB composite generation function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f592790b-26cf-465c-9511-6c751799f39f", + "metadata": {}, + "outputs": [], + "source": [ + "cutouts = [cutout[band] for band in bands]\n", + "multibandexposure = MultibandExposure.fromExposures(bands, cutouts)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "23848e70-3f5d-40e4-bc6b-0e208af3f044", + "metadata": {}, + "outputs": [], + "source": [ + "showRGB(multibandexposure.image, bgr='gri', figsize=(3,3), stretch=60, Q=10, name=name)" + ] + }, + { + "cell_type": "markdown", + "id": "38a9221a-1ecf-48f7-9584-991f9d8b512c", + "metadata": {}, + "source": [ + "We can use a more recent method for making color composite images, that was used in the First Look release! Here's an alternative RGB function:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bb37f19f-cb5f-4fb0-b2ef-7ee2308104ba", + "metadata": {}, + "outputs": [], + "source": [ + "from lsst.pipe.tasks.prettyPictureMaker import PrettyPictureConfig, PrettyPictureTask\n", + "from lsst.pipe.tasks.prettyPictureMaker._task import ChannelRGBConfig\n", + "\n", + "def prettyRGB(images, ax=None, fp=None, figsize=(8,8), stretch=250, Q=0.7, name=None):\n", + " \"\"\"Display an RGB (irg) color composite image with matplotlib using the prettyPictureMaker pipe task.\n", + " \n", + " Parameters\n", + " ----------\n", + " images : dict\n", + " Dictionary of images to display.\n", + " ax : `matplotlib.axes.Axes`\n", + " Axis in a `matplotlib.Figure` to display the image.\n", + " If `axis` is `None` then a new figure is created.\n", + " figsize: tuple\n", + " Size of the `matplotlib.Figure` created.\n", + " If `ax` is not `None` then this parameter is ignored.\n", + " stretch: int\n", + " The linear stretch of the image.\n", + " Q: int\n", + " The Asinh softening parameter.\n", + " name: str\n", + " The name of the object/field to be displayed.\n", + " \"\"\"\n", + "\n", + " prettyPicConfig = PrettyPictureTask.ConfigClass()\n", + " # Magic from Nate Lust:\n", + " prettyPicConfig.localContrastConfig.doLocalContrast = False\n", + " prettyPicConfig.localContrastConfig.sigma = 30\n", + " prettyPicConfig.localContrastConfig.clarity = 0.8\n", + " prettyPicConfig.localContrastConfig.shadows = 0\n", + " prettyPicConfig.localContrastConfig.highlights = -1.5\n", + " prettyPicConfig.localContrastConfig.maxLevel = 2\n", + " prettyPicConfig.imageRemappingConfig.absMax = 11000\n", + " prettyPicConfig.luminanceConfig.max = 100\n", + " prettyPicConfig.luminanceConfig.stretch = stretch # from kwargs\n", + " prettyPicConfig.luminanceConfig.floor = 0\n", + " prettyPicConfig.luminanceConfig.Q = Q # from kwargs\n", + " prettyPicConfig.luminanceConfig.highlight = 0.905882\n", + " prettyPicConfig.luminanceConfig.shadow = 0.12\n", + " prettyPicConfig.luminanceConfig.midtone = 0.25\n", + " prettyPicConfig.doPSFDeconcovlve = False # sic\n", + " prettyPicConfig.exposureBrackets = None\n", + " prettyPicConfig.colorConfig.maxChroma = 80\n", + " prettyPicConfig.colorConfig.saturation = 0.6\n", + " prettyPicConfig.cieWhitePoint = (0.28, 0.28)\n", + " prettyPicConfig.channelConfig = dict(\n", + " g=ChannelRGBConfig(r=0.0, g=0.0, b=1.0),\n", + " r=ChannelRGBConfig(r=0.0, g=1.0, b=0.0),\n", + " i=ChannelRGBConfig(r=1.0, g=0.0, b=0.0),\n", + " )\n", + " prettyPicTask = PrettyPictureTask(config=prettyPicConfig)\n", + " \n", + " bands = \"gri\"\n", + " coaddG = images['g']\n", + " coaddR = images['r']\n", + " coaddI = images['i']\n", + " \n", + " prettyPicInputs = prettyPicTask.makeInputsFromExposures(i=coaddI, r=coaddR, g=coaddG)\n", + " coaddRgbStruct = prettyPicTask.run(prettyPicInputs)\n", + " coaddRgb = coaddRgbStruct.outputRGB\n", + "\n", + " if ax is None:\n", + " fig = plt.figure(figsize=figsize)\n", + " ax = fig.add_subplot(1,1,1)\n", + " \n", + " plt.axis(\"off\")\n", + " ax.imshow(coaddRgb, interpolation='nearest', origin='lower')\n", + "\n", + " if name is not None:\n", + " plt.text(0, 31, name, color='white', fontsize=12, horizontalalignment='left', verticalalignment='top')\n", + " \n", + " plt.text(0, 2, 'PrettyPictureTask '+bands, color='white', fontsize=12, horizontalalignment='left', verticalalignment='top')\n", + " \n", + " plt.tight_layout();\n", + "\n", + " return" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dc1f388a-18d9-473d-bd5c-994863f324f5", + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(8,4))\n", + "\n", + "ax = fig.add_subplot(1,2,1)\n", + "showRGB(multibandexposure.image, bgr='gri', ax=ax, stretch=60, Q=10, name=name)\n", + "\n", + "ax = fig.add_subplot(1,2,2)\n", + "prettyRGB(cutout, ax=ax, stretch=750, Q=0.7, name=name)" + ] + }, + { + "cell_type": "markdown", + "id": "37c0f387-d8a9-4723-9668-f227703fa105", + "metadata": {}, + "source": [ + "The `prettyPictureMaker` code does a better job at bringing out the color contrast.\n", + "\n", + "Choosing the stretch and Q, in either method, can be a bit fiddly - this is best done when visualizing the whole set of cutouts in a gallery. This is what we will do next." + ] + }, + { + "cell_type": "markdown", + "id": "173ea5ab-2f3c-4a46-a65a-84430244f966", + "metadata": {}, + "source": [ + "## Multiple Objects, Multiple Bands, Gallery of Images\n", + "\n", + "What we want is a function that takes in a table of targets, makes cutouts for them all (if they haven't been made already), and then enables visualization of them as a gallery of RGB color composite images (all on the same stretch)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ecefadb1-5401-4573-a91f-3861b8f793bc", + "metadata": {}, + "outputs": [], + "source": [ + "class StampCollector():\n", + " def __init__(self):\n", + " self.bands = ['u','g','r','i','z','y']\n", + " self.cutoutSize = geom.ExtentI(32, 32)\n", + " self.cutouts = {}\n", + " butler = Butler(\"dp1\", collections=\"LSSTComCam/DP1\")\n", + " assert butler is not None\n", + " butler.get_dataset_type('deep_coadd')\n", + " self.skymap = butler.get(\"skyMap\")\n", + " return\n", + " \n", + " def get_cutouts(self, targets, bands=None):\n", + " \"\"\"Read in a table of target RA and Dec positions (in degrees) and make cutout images for them all\n", + " \n", + " Parameters\n", + " ----------\n", + " targets : `astropy.Table`\n", + " Table of target coordinates.\n", + " bands : list of strings\n", + " Which bands to make cutouts in.\n", + " \"\"\"\n", + " self.targets = targets\n", + " # Make IAU names if they don't exist already:\n", + " if 'name' not in self.targets.columns:\n", + " print(\"Adding IAU names...\")\n", + " self.add_iau_names()\n", + "\n", + " if bands is not None:\n", + " self.bands = bands\n", + "\n", + " print(\"Making cutouts:\")\n", + "\n", + " # Loop over targets and bands, extracting cutouts:\n", + " for i in range(len(self.targets)):\n", + "\n", + " print(\" \", self.targets['name'][i],\": \",end=\"\")\n", + "\n", + " radec = geom.SpherePoint(self.targets['ra'][i], self.targets['dec'][i], geom.degrees) \n", + " tractInfo = self.skymap.findTract(radec)\n", + " patchInfo = tractInfo.findPatch(radec)\n", + " patch = tractInfo.getSequentialPatchIndex(patchInfo)\n", + " tract = tractInfo.getId()\n", + " \n", + " xy = geom.PointI(tractInfo.getWcs().skyToPixel(radec))\n", + " bbox = geom.BoxI(xy - cutoutSize // 2, cutoutSize)\n", + " parameters = {'bbox': bbox}\n", + "\n", + " if not self.targets['name'][i] in self.cutouts:\n", + " self.cutouts[self.targets['name'][i]] = {}\n", + " \n", + " for band in self.bands:\n", + "\n", + " if not band in self.cutouts[self.targets['name'][i]]:\n", + " dataId = {'tract': tract, 'patch': patch, 'band': band}\n", + " try:\n", + " self.cutouts[self.targets['name'][i]][band] = butler.get(\"deep_coadd\", parameters=parameters, dataId=dataId)\n", + " print(band,end=\"\")\n", + " except:\n", + " print(\".\",end=\"\")\n", + " else:\n", + " print(\".\",end=\"\")\n", + "\n", + " print()\n", + " \n", + " return\n", + " \n", + " def add_iau_names(self):\n", + " names = []\n", + " for i in range(len(self.targets)):\n", + " rahms = self.targets['ra'][i] * u.hourangle\n", + " decdms = self.targets['dec'][i] * u.deg\n", + " coordinates = SkyCoord(ra=rahms, dec=decdms, frame='icrs')\n", + " names.append((f'EUCLID J{coordinates.ra.to_string(unit=u.hourangle, sep=\"\", precision=1, pad=True)}'\n", + " f'{coordinates.dec.to_string(sep=\"\", precision=0, alwayssign=True, pad=True)}'))\n", + " self.targets.add_column(names, name='name', index=0)\n", + " return\n", + " \n", + " def make_gallery(self, nx=3, style='Pretty', bgr=\"gri\", stretch=750, Q=0.7):\n", + " \"\"\"Display a gallery of RGB color composite images with matplotlib.\n", + " \n", + " Parameters\n", + " ----------\n", + " nx : integer\n", + " Number of images in a row of the gallery.\n", + " style : str\n", + " Method to use, Lupton (astropy) or Pretty (Rubin)\n", + " bgr : sequence\n", + " A 3-element sequence of filter names (i.e. keys of the exps dict) indicating what band\n", + " to use for each channel. If `image` only has three filters then this parameter is ignored\n", + " and the filters in the image are used.\n", + " stretch: int\n", + " The linear stretch of the image.\n", + " Q: int\n", + " The Asinh softening parameter.\n", + " \"\"\"\n", + " # Figure out gallery dimensions:\n", + " rem = len(self.targets) % nx\n", + " ny = int((len(self.targets) - rem) / nx) + 1\n", + " print(\"Making gallery, \",nx,\" wide by \",ny,\" deep.\")\n", + "\n", + " # Create an nx by ny grid of subplots\n", + " width = 12\n", + " fig = plt.figure(figsize=(width, ny*(width/nx)))\n", + " plt.axis(\"off\")\n", + "\n", + " # Loop over targets, stepping through subplots:\n", + " for k in range(len(self.targets)):\n", + " ax = fig.add_subplot(ny,nx,k+1)\n", + " if style == 'Lupton':\n", + " self.show_RGB(self.targets['name'][k], bgr=bgr, ax=ax, stretch=stretch, Q=Q)\n", + " else:\n", + " assert bgr == \"gri\"\n", + " self.pretty_RGB(self.targets['name'][k], ax=ax, stretch=stretch, Q=Q)\n", + " \n", + " plt.tight_layout()\n", + " \n", + " return\n", + "\n", + " def show_RGB(self, name, bgr=\"gri\", ax=None, fp=None, figsize=(8,8), stretch=100, Q=5):\n", + " \"\"\"Display an RGB color composite image with matplotlib.\n", + " \n", + " Parameters\n", + " ----------\n", + " name : str\n", + " Name of the object/field to display.\n", + " bgr : sequence\n", + " A 3-element sequence of filter names (i.e. keys of the exps dict) indicating what band\n", + " to use for each channel. If `image` only has three filters then this parameter is ignored\n", + " and the filters in the image are used.\n", + " ax : `matplotlib.axes.Axes`\n", + " Axis in a `matplotlib.Figure` to display the image.\n", + " If `axis` is `None` then a new figure is created.\n", + " figsize: tuple\n", + " Size of the `matplotlib.Figure` created.\n", + " If `ax` is not `None` then this parameter is ignored.\n", + " stretch: int\n", + " The linear stretch of the image.\n", + " Q: int\n", + " The Asinh softening parameter.\n", + " \"\"\"\n", + " bands = [bgr[0],bgr[1],bgr[2]]\n", + " cutouts = [self.cutouts[name][band] for band in bands]\n", + " image = MultibandExposure.fromExposures(bands, cutouts).image\n", + " \n", + " # Extract the primary image component of each Exposure with the .image property, and use .array to get a NumPy array view.\n", + " rgb = make_lupton_rgb(image_r=image[bgr[2]].array, # numpy array for the r channel\n", + " image_g=image[bgr[1]].array, # numpy array for the g channel\n", + " image_b=image[bgr[0]].array, # numpy array for the b channel\n", + " stretch=stretch, Q=Q) # parameters used to stretch and scale the pixel values\n", + " if ax is None:\n", + " fig = plt.figure(figsize=figsize)\n", + " ax = fig.add_subplot(1,1,1)\n", + " \n", + " plt.axis(\"off\")\n", + " ax.imshow(rgb, interpolation='nearest', origin='lower')\n", + " \n", + " if name is not None:\n", + " plt.text(0, 31, name, color='white', fontsize=16, horizontalalignment='left', verticalalignment='top')\n", + " \n", + " plt.text(0, 2, bgr, color='white', fontsize=16, horizontalalignment='left', verticalalignment='top')\n", + "\n", + " del image, cutouts, rgb\n", + " return\n", + "\n", + " def pretty_RGB(self, name, ax=None, fp=None, figsize=(8,8), stretch=750, Q=0.7):\n", + " \"\"\"Display an RGB (irg) color composite image with matplotlib using the prettyPictureMaker pipe task.\n", + " \n", + " Parameters\n", + " ----------\n", + " name : str\n", + " Name of the object/field to display.\n", + " ax : `matplotlib.axes.Axes`\n", + " Axis in a `matplotlib.Figure` to display the image.\n", + " If `axis` is `None` then a new figure is created.\n", + " figsize: tuple\n", + " Size of the `matplotlib.Figure` created.\n", + " If `ax` is not `None` then this parameter is ignored.\n", + " stretch: int\n", + " The linear stretch of the image.\n", + " Q: int\n", + " The Asinh softening parameter.\n", + " name: str\n", + " The name of the object/field to be displayed.\n", + " \"\"\"\n", + " \n", + " prettyPicConfig = PrettyPictureTask.ConfigClass()\n", + " # Magic from Nate Lust:\n", + " prettyPicConfig.localContrastConfig.doLocalContrast = False\n", + " prettyPicConfig.localContrastConfig.sigma = 30\n", + " prettyPicConfig.localContrastConfig.clarity = 0.8\n", + " prettyPicConfig.localContrastConfig.shadows = 0\n", + " prettyPicConfig.localContrastConfig.highlights = -1.5\n", + " prettyPicConfig.localContrastConfig.maxLevel = 2\n", + " prettyPicConfig.imageRemappingConfig.absMax = 11000\n", + " prettyPicConfig.luminanceConfig.max = 100\n", + " prettyPicConfig.luminanceConfig.stretch = stretch # from kwargs\n", + " prettyPicConfig.luminanceConfig.floor = 0\n", + " prettyPicConfig.luminanceConfig.Q = Q # from kwargs\n", + " prettyPicConfig.luminanceConfig.highlight = 0.905882\n", + " prettyPicConfig.luminanceConfig.shadow = 0.12\n", + " prettyPicConfig.luminanceConfig.midtone = 0.25\n", + " prettyPicConfig.doPSFDeconcovlve = False # sic\n", + " prettyPicConfig.exposureBrackets = None\n", + " prettyPicConfig.colorConfig.maxChroma = 80\n", + " prettyPicConfig.colorConfig.saturation = 0.6\n", + " prettyPicConfig.cieWhitePoint = (0.28, 0.28)\n", + " prettyPicConfig.channelConfig = dict(\n", + " g=ChannelRGBConfig(r=0.0, g=0.0, b=1.0),\n", + " r=ChannelRGBConfig(r=0.0, g=1.0, b=0.0),\n", + " i=ChannelRGBConfig(r=1.0, g=0.0, b=0.0),\n", + " )\n", + " prettyPicTask = PrettyPictureTask(config=prettyPicConfig)\n", + " \n", + " bands = \"gri\"\n", + " coaddG = self.cutouts[name]['g']\n", + " coaddR = self.cutouts[name]['r']\n", + " coaddI = self.cutouts[name]['i']\n", + " \n", + " prettyPicInputs = prettyPicTask.makeInputsFromExposures(i=coaddI, r=coaddR, g=coaddG)\n", + " coaddRgbStruct = prettyPicTask.run(prettyPicInputs)\n", + " coaddRgb = coaddRgbStruct.outputRGB\n", + " \n", + " if ax is None:\n", + " fig = plt.figure(figsize=figsize)\n", + " ax = fig.add_subplot(1,1,1)\n", + " \n", + " plt.axis(\"off\")\n", + " ax.imshow(coaddRgb, interpolation='nearest', origin='lower')\n", + " \n", + " if name is not None:\n", + " plt.text(0, 31, name, color='white', fontsize=12, horizontalalignment='left', verticalalignment='top')\n", + " \n", + " plt.text(0, 2, bands, color='white', fontsize=12, horizontalalignment='left', verticalalignment='top')\n", + " \n", + " plt.tight_layout();\n", + " \n", + " return" + ] + }, + { + "cell_type": "markdown", + "id": "3c17463f-5e08-4cd9-a5d4-4a832cb2b430", + "metadata": {}, + "source": [ + "With our `StampCollector` ready, we can prepare and provide a table of targets. Here's a simple example with just two Euclid lens candidates. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "47fdc1aa-68ca-486d-b240-a2a090dfd702", + "metadata": {}, + "outputs": [], + "source": [ + "target_table = Table(names=('ra', 'dec'), dtype=('f4', 'f4'))\n", + "target_table.add_row((59.496380, -48.494726))\n", + "target_table.add_row((59.626134, -49.06175))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3f96ccfd-4018-43bd-81bb-0b757d5c8475", + "metadata": {}, + "outputs": [], + "source": [ + "stamp_collector = StampCollector()\n", + "stamp_collector.get_cutouts(target_table,bands=['g','r','i'])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b5517c76-d685-4d51-b455-4061a6e5672a", + "metadata": {}, + "outputs": [], + "source": [ + "stamp_collector.targets" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b6424825-9bde-4db7-89e4-1a55f307ad4b", + "metadata": {}, + "outputs": [], + "source": [ + "stamp_collector.make_gallery(stretch=750,Q=0.7)" + ] + }, + { + "cell_type": "markdown", + "id": "dfee8c96-4229-4663-bcb7-b25bd3ab3491", + "metadata": {}, + "source": [ + "Now let's get all the Euclid lens candidates, in all DP1 fields. These can be downloaded as CSV from a Google sheet on the web, using the pandas native `read_csv` function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "356f0bba-95f5-4755-aa90-5f0e10e8e9ae", + "metadata": {}, + "outputs": [], + "source": [ + "import pandas as pd\n", + "\n", + "sheet_id = '1OInIecou_c2NqVdpTBSYB2V_j2vanKvKrTBf0qtLCpE' \n", + "csv_url = f'https://docs.google.com/spreadsheets/d/{sheet_id}/export?format=csv' \n", + "\n", + "# Read the data into a pandas DataFrame\n", + "df = pd.read_csv(csv_url)\n", + "\n", + "# Convert the DataFrame to an Astropy Table\n", + "all_targets = Table.from_pandas(df)[['Euclid_RA', 'Euclid_Dec']] \n", + "all_targets.rename_column('Euclid_RA', 'ra')\n", + "all_targets.rename_column('Euclid_Dec', 'dec')\n", + "print(all_targets)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f07a7260-737f-4738-8b71-f49fcfe90a3c", + "metadata": {}, + "outputs": [], + "source": [ + "all_stamp_collector = StampCollector()\n", + "all_stamp_collector.get_cutouts(all_targets, bands=['g','r','i'])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c4dd572c-53ce-4ed7-99ff-daa75a6c65c6", + "metadata": {}, + "outputs": [], + "source": [ + "all_stamp_collector.make_gallery(nx=3,style=\"Pretty\",stretch=750,Q=0.7)" + ] + }, + { + "cell_type": "markdown", + "id": "c84d8a66-f711-4bac-96bc-ba3b8e0579c6", + "metadata": {}, + "source": [ + "## Further Work\n", + "\n", + "* Investigate the Cutout Factory approach - if we can get the world coordinate to image coordinate transformation right, this should give us a bit of a speed up for cases where we have multiple targets in the same patch.\n", + "* The gallery needs work: better RGB representation (try and follow the First Look set up?), and a more efficient packing of small images (including scaled fontsize for the overlays).\n", + "* Add a method to package ugrizy cutouts into the right format for the Rubin SharPy - and then run that code to sharpen up the images and look for lensed arcs!" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "LSST", + "language": "python", + "name": "lsst" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}