diff --git a/DP1/100_How_to_Use_RSP_Tools/103_Image_access_and_display/103_4_Small_Image_Cutout.ipynb b/DP1/100_How_to_Use_RSP_Tools/103_Image_access_and_display/103_4_Small_Image_Cutout.ipynb index cd98dae..aece252 100644 --- a/DP1/100_How_to_Use_RSP_Tools/103_Image_access_and_display/103_4_Small_Image_Cutout.ipynb +++ b/DP1/100_How_to_Use_RSP_Tools/103_Image_access_and_display/103_4_Small_Image_Cutout.ipynb @@ -22,7 +22,7 @@ "Data Release: DP1
\n", "Container Size: large
\n", "LSST Science Pipelines version: r29.2.0
\n", - "Last verified to run: 2025-12-16
\n", + "Last verified to run: 2026-01-15
\n", "Repository: github.com/lsst/tutorial-notebooks
" ] }, @@ -100,11 +100,9 @@ "import matplotlib.pyplot as plt\n", "\n", "import lsst.afw.display as afwDisplay\n", - "import lsst.afw.image as afwImage\n", - "from lsst.afw.image import ExposureF\n", + "from lsst.afw.image import ExposureF, ImageF, MaskedImageF\n", "from lsst.afw.fits import MemFileManager\n", "\n", - "\n", "from lsst.rsp.utils import get_pyvo_auth\n", "from lsst.rsp.service import get_siav2_service\n", "\n", @@ -115,7 +113,7 @@ "from astropy import units as u\n", "from astropy.coordinates import Angle\n", "from astropy.time import Time\n", - "\n", + "from astropy.wcs import WCS\n", "from astropy.io import fits\n", "import io\n" ] @@ -328,9 +326,9 @@ "id": "d945e191-7ee4-458a-bb0f-fb00cd219e30", "metadata": {}, "source": [ - "Lastly, call the Rubin Image Cutout Service. This notebook demonstrates two types of cutout services: `cutout-sync` which, by default, returns just the image extension of the LSST imaging; `cutout-sync-exposure`, which returns the full set of metadata and image extensions that are contained in the `ExposureF` data type. This section will demonstrate the use of `cutout-sync-exposure` and Section 4 below will demonstrate the other. \n", + "Lastly, call the Rubin Image Cutout Service. This notebook demonstrates three types of cutout services: `cutout-sync` which, by default, returns just the image extension of the LSST imaging; `cutout-sync-exposure`, which returns the full set of metadata and image extensions that are contained in the `ExposureF` data type, and `cutout-sync-mask` which returns the mask extension with the science and variance images (`MaskedImageF` data type). This section will demonstrate the use of the default cutout service, `cutout-sync`, and Section 4 below will demonstrate the other two. \n", "\n", - "To use the cutout service in this example, the IVOA procedure `cutout-sync-exposure` is called using `get_adhocservice_by_id`. It is done by feeding the data link created above (called `dl_result`) to `from_resource`. Since the Rubin DP1 imaging is proprietary it is necessary to again provide the authorization for the current RSP session. Do this using the `get_pyvo_auth` function." + "To use the cutout service in this example, the IVOA procedure `cutout-sync` is called using `get_adhocservice_by_id`. It is done by feeding the data link created above (called `dl_result`) to `from_resource`. Since the Rubin DP1 imaging is proprietary it is necessary to again provide the authorization for the current RSP session. Do this using the `get_pyvo_auth` function." ] }, { @@ -341,7 +339,7 @@ "outputs": [], "source": [ "sq = SodaQuery.from_resource(dl_result,\n", - " dl_result.get_adhocservice_by_id(\"cutout-sync-exposure\"),\n", + " dl_result.get_adhocservice_by_id(\"cutout-sync\"),\n", " session=get_pyvo_auth())" ] }, @@ -354,7 +352,9 @@ "\n", "### 3.1. Define cutout center and edge\n", "\n", - "Only two shape definitions are supported: a circle function, and a polygon function can be used to define the cutout dimensions. These shape definitions do not produce circle or polygon cutouts, but rather are methods for defining the edges of cutouts with 4 sides. In the case of circle, the resulting cutout is always a square, with edge size that is the same as the circle diameter. In the case of a polygon, either a square or a rectangular cutout will result, depending on whether the length and width edge dimensions are different values. Only cutouts with 4 corners and 90 degree angles are supported. " + "Only two shape definitions are supported: a circle function, and a polygon function can be used to define the cutout dimensions. These shape definitions do not produce circle or polygon cutouts, but rather are methods for defining the edges of cutouts with 4 sides. In the case of circle, the resulting cutout is always a square, with edge size that is the same as the circle diameter. In the case of a polygon, either a square or a rectangular cutout will result, depending on whether the length and width edge dimensions are different values. Only cutouts with 4 corners and 90 degree angles are supported. \n", + "\n", + "`.circle` defaults to assuming the units are degrees; this notebook demonstrates its use when specifying the units with astropy." ] }, { @@ -364,21 +364,25 @@ "metadata": {}, "outputs": [], "source": [ - "spherePoint = geom.SpherePoint(target_ra*geom.degrees, target_dec*geom.degrees)\n", + "cutout_ra = target_ra * u.deg\n", + "cutout_dec = target_dec * u.deg\n", + "\n", "Radius = 0.01 * u.deg\n", - "sq.circle = (spherePoint.getRa().asDegrees() * u.deg,\n", - " spherePoint.getDec().asDegrees() * u.deg,\n", - " Radius)" + "\n", + "sq.circle = (cutout_ra, cutout_dec, Radius)\n", + "\n", + "cutout_bytes = sq.execute_stream().read()\n", + "sq.raise_if_error()" ] }, { "cell_type": "markdown", - "id": "30640af4-a75c-4132-8f2e-692fcf9694eb", + "id": "9be31876-1afa-4610-aaec-cc3497ad6930", "metadata": {}, "source": [ - "### 3.2. Retrieve the cutout\n", + "### 3.2. Retrieve the cutout with LSST\n", "\n", - "The recommended method is to recieve the cutout into memory in the `ExposureF` format." + "The cutout can be read into memory in the LSST pipelines `ImageF` format, but does not include the WCS information." ] }, { @@ -388,11 +392,9 @@ "metadata": {}, "outputs": [], "source": [ - "cutout_bytes = sq.execute_stream().read()\n", - "sq.raise_if_error()\n", "mem = MemFileManager(len(cutout_bytes))\n", "mem.setData(cutout_bytes, len(cutout_bytes))\n", - "exposure = ExposureF(mem)\n" + "image = ImageF(mem)" ] }, { @@ -400,7 +402,7 @@ "id": "324be87a-337b-4db7-b57d-842f14614613", "metadata": {}, "source": [ - "Display the image extension of the `ExposureF` cutout." + "Display the `ImageF` cutout." ] }, { @@ -412,7 +414,7 @@ "source": [ "display = afwDisplay.Display()\n", "display.scale('asinh', 'zscale')\n", - "display.image(exposure.image)\n", + "display.image(image)\n", "plt.show()" ] }, @@ -421,7 +423,7 @@ "id": "458a4d25-ca95-4f92-b76e-986b7a03b320", "metadata": {}, "source": [ - "> **Figure 1:** The cutout image, displayed in grayscale with a scale bar at right.\n", + "> **Figure 1:** The cutout image, displayed in pixel coordinates using LSST pipeline tools in grayscale with a scale bar at right.\n", "\n", "#### 3.2.1. Option to save cutout to disk\n", "\n", @@ -429,7 +431,7 @@ "\n", "The cutouts saved to disk with the following commands will be stored in a temporary folder in the user's home directory.\n", "\n", - "After saving the cutout as a FITS file, read it in as an `ExposureF` and display it." + "After saving the cutout as a FITS file, read it in as an `ImageF` and display it." ] }, { @@ -450,10 +452,10 @@ "# with open(sodaCutout, 'bw') as f:\n", "# f.write(sq.execute_stream().read())\n", "\n", - "# cutout = ExposureF(sodaCutout)\n", + "# cutout = ImageF(sodaCutout)\n", "# display = afwDisplay.Display()\n", "# display.scale('asinh', 'zscale')\n", - "# display.image(cutout.image)\n", + "# display.image(cutout)\n", "# plt.show()" ] }, @@ -463,23 +465,96 @@ "metadata": {}, "source": [ "Future planned options for the Rubin cutout service, including the potential to retrieve other image formats such as jpeg, are listed at the Rubin Science Platform image cutout implementation strategy document.\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "30640af4-a75c-4132-8f2e-692fcf9694eb", + "metadata": {}, + "source": [ + "### 3.3. Retrieve the cutout with astropy\n", "\n", - "**Note:** `ExposureF` contains multiple Header-Data Units (HDUs) with extensive exposure metadata. To reduce file size, one option is to save only the image, mask, and variance planes, along with the WCS recorded in the Primary HDU. This can be done by writing `exposure.maskedImage`-a three-plane object containing the image, mask, and variance-into a FITS file using `writeFitsImage`." + "The format of the default cutout is a 2-extension fits: the primary (extension 0) header data unit (HDU) has the header describing the data set, and the secondary HDU (extension 1) has only the header that describes the pixel data (i.e. science image).\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "de03af16-26d1-4268-998b-43f33bb35a67", + "id": "6cd941c5-7e03-4c1b-a327-03f5ebc578f1", "metadata": {}, "outputs": [], "source": [ - "# fits.info(sodaCutout)\n", + "hdul = fits.open(io.BytesIO(cutout_bytes))\n" + ] + }, + { + "cell_type": "markdown", + "id": "87c5e852-f42c-403f-b357-130a69b2b2e5", + "metadata": {}, + "source": [ + "Explore the format of the returned HDU extensions. There are two HDUs, where the second one labeled Image (extension 1) contains the image pixels and WCS header." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "00b1efae-7b59-422b-9708-1036120c4d8d", + "metadata": {}, + "outputs": [], + "source": [ + "hdul.info()\n", + "primary_hdu = hdul[0]\n", + "image_hdu = hdul[1]\n" + ] + }, + { + "cell_type": "markdown", + "id": "a5a6be22-ebe9-4fab-ac3a-cecd2126994b", + "metadata": {}, + "source": [ + "Below, store the headers to allow the option to print (note that the primary header is very long)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "da4cc3eb-16e3-44c7-9cf7-0d8214bc0671", + "metadata": {}, + "outputs": [], + "source": [ + "primary_header = hdul[0].header\n", + "image_header = hdul[1].header\n" + ] + }, + { + "cell_type": "markdown", + "id": "88be8f44-3226-4fda-97eb-ea0b272a88ec", + "metadata": {}, + "source": [ + "Plot the cutout using astropy tools, with WCS info from the image header used to plot in sky coordinates." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0b755176-ce5f-4238-9cd9-4798a0dfe133", + "metadata": {}, + "outputs": [], + "source": [ + "wcs = WCS(image_header)\n", "\n", - "# sodaCutout_small = os.path.join(tempdir, 'cutout-circle_small.fits')\n", - "# afwDisplay.writeFitsImage(sodaCutout_small,\n", - "# exposure.maskedImage,\n", - "# wcs=exposure.wcs)" + "fig = plt.figure(figsize=(6, 6))\n", + "ax = fig.add_subplot(111, projection=wcs)\n", + "ax.imshow(image_hdu.data, cmap='gray', norm='asinh', origin='lower')\n" + ] + }, + { + "cell_type": "markdown", + "id": "ca5f197f-24a9-4440-8c8e-9f3398ebea26", + "metadata": {}, + "source": [ + "> **Figure 2:** The cutout image, displayed in sky coordinates using astropy.\n" ] }, { @@ -487,13 +562,13 @@ "id": "a1af5867-6da1-40f0-99ce-703938338ce4", "metadata": {}, "source": [ - "### 3.3. Use polygon to define shape\n", + "### 3.4. Use polygon to define shape\n", "\n", "It is also possible to define the cutout geometry using a polygon, which enables the cutout to be rectangular, but not necessarily be square. For this, use `polygon`, which takes as input the four corners in celestial coordinates. A minimum of 3 vertices are required (the line from the last vertex back to the first is implicit) Vertices must be ordered in the counter-clockwise direction. For example: a polygon is defined as a set of 4 (x,y) coordinates from (12,34) to (14,34) to (14,36) to (12,36) and (implicitly) back to (12,34) as:\n", "\n", "POLYGON=12 34 14 34 14 36 12 36\n", "\n", - "Since the center of the cutout is already defined in RA and Dec in the cells above (`spherePoint`), this example will define each x,y set as RA+/-edge and Dec+/-edge.\n", + "Since the center of the cutout is already defined in unit degrees in the cells above, this example will define each x,y set as RA+/-edge and Dec+/-edge.\n", "\n", "> **Warning:** Visit images are rotated, and although it is the \"Declination edge\" that is defined to be smaller, that corresponds to the x-axis when the cutout is displayed below, due to image rotation." ] @@ -506,26 +581,26 @@ "outputs": [], "source": [ "sqp = SodaQuery.from_resource(dl_result,\n", - " dl_result.get_adhocservice_by_id(\"cutout-sync-exposure\"),\n", + " dl_result.get_adhocservice_by_id(\"cutout-sync\"),\n", " session=get_pyvo_auth())\n", "\n", "ra_edge = 0.02 * u.deg\n", "de_edge = 0.005 * u.deg\n", "\n", - "sqp.polygon = (spherePoint.getRa().asDegrees() * u.deg - ra_edge,\n", - " spherePoint.getDec().asDegrees() * u.deg - de_edge,\n", - " spherePoint.getRa().asDegrees() * u.deg - ra_edge,\n", - " spherePoint.getDec().asDegrees() * u.deg + de_edge,\n", - " spherePoint.getRa().asDegrees() * u.deg + ra_edge,\n", - " spherePoint.getDec().asDegrees() * u.deg + de_edge,\n", - " spherePoint.getRa().asDegrees() * u.deg + ra_edge,\n", - " spherePoint.getDec().asDegrees() * u.deg - de_edge)\n", + "sqp.polygon = (cutout_ra - ra_edge,\n", + " cutout_dec - de_edge,\n", + " cutout_ra - ra_edge,\n", + " cutout_dec + de_edge,\n", + " cutout_ra + ra_edge,\n", + " cutout_dec + de_edge,\n", + " cutout_ra + ra_edge,\n", + " cutout_dec - de_edge)\n", "\n", "cutout_bytes = sqp.execute_stream().read()\n", "sqp.raise_if_error()\n", "mem = MemFileManager(len(cutout_bytes))\n", "mem.setData(cutout_bytes, len(cutout_bytes))\n", - "polygon = ExposureF(mem)" + "polygon = ImageF(mem)" ] }, { @@ -537,7 +612,7 @@ "source": [ "display = afwDisplay.Display()\n", "display.scale('asinh', 'zscale')\n", - "display.mtv(polygon.image)\n", + "display.mtv(polygon)\n", "plt.show()" ] }, @@ -546,7 +621,7 @@ "id": "c83e5413-2edf-41fd-8542-b622ba8324b0", "metadata": {}, "source": [ - "> **Figure 2:** A rectangular polygon cutout from a rotated `visit_image`." + "> **Figure 3:** A rectangular polygon cutout from a rotated `visit_image`." ] }, { @@ -554,7 +629,7 @@ "id": "9041e6f0-e9bc-46df-a9b8-ad0c61deb6bb", "metadata": {}, "source": [ - "### 3.4. Correcting for cos(d)\n", + "### 3.5. Correcting for cos(d)\n", "\n", "There is an important difference to note between the circle and polygon shape definitions. The angular distance on the sky that defines the circular cutout size already accounts for the difference in angular distance in the RA direction is smaller by a factor of cos(declination), where declination is in units radians. The difference increases with higher declination. However, the polygon definition does not automatically account for this cosine factor. Thus, circle and polygon cutout definitions using the same cutout edge length will not match size in the RA direction. The 2 cells below demonstrate how to make this correction to the polygon cutout definition to create symmetric cutouts with polygon. Here, reset the edge sizes to be the same as `Radius` from the circle definition above.\n", "\n", @@ -569,24 +644,21 @@ "outputs": [], "source": [ "sq2 = SodaQuery.from_resource(dl_result,\n", - " dl_result.get_adhocservice_by_id(\"cutout-sync-exposure\"),\n", + " dl_result.get_adhocservice_by_id(\"cutout-sync\"),\n", " session=get_pyvo_auth())\n", "\n", - "edge1 = Radius\n", - "edge2 = Radius\n", - "\n", - "sq2.polygon = (spherePoint.getRa().asDegrees() * u.deg - edge1,\n", - " spherePoint.getDec().asDegrees() * u.deg - edge2,\n", - " spherePoint.getRa().asDegrees() * u.deg - edge1,\n", - " spherePoint.getDec().asDegrees() * u.deg + edge2,\n", - " spherePoint.getRa().asDegrees() * u.deg + edge1,\n", - " spherePoint.getDec().asDegrees() * u.deg + edge2)\n", + "sq2.polygon = (cutout_ra - Radius,\n", + " cutout_dec - Radius,\n", + " cutout_ra - Radius,\n", + " cutout_dec + Radius,\n", + " cutout_ra + Radius,\n", + " cutout_dec + Radius)\n", "\n", "cutout_bytes = sq2.execute_stream().read()\n", "sq2.raise_if_error()\n", "mem = MemFileManager(len(cutout_bytes))\n", "mem.setData(cutout_bytes, len(cutout_bytes))\n", - "exposure2 = ExposureF(mem)" + "polygon2 = ImageF(mem)" ] }, { @@ -604,27 +676,28 @@ "metadata": {}, "outputs": [], "source": [ + "spherePoint = geom.SpherePoint(target_ra*geom.degrees, target_dec*geom.degrees)\n", "a = Angle(spherePoint.getDec().asDegrees(), u.deg)\n", "cosd = np.cos(a.radian)\n", "\n", "sq3 = SodaQuery.from_resource(dl_result,\n", - " dl_result.get_adhocservice_by_id(\"cutout-sync-exposure\"),\n", + " dl_result.get_adhocservice_by_id(\"cutout-sync\"),\n", " session=get_pyvo_auth())\n", "\n", - "sq3.polygon = (spherePoint.getRa().asDegrees() * u.deg - edge1/cosd,\n", - " spherePoint.getDec().asDegrees() * u.deg - edge2,\n", - " spherePoint.getRa().asDegrees() * u.deg - edge1/cosd,\n", - " spherePoint.getDec().asDegrees() * u.deg + edge2,\n", - " spherePoint.getRa().asDegrees() * u.deg + edge1/cosd,\n", - " spherePoint.getDec().asDegrees() * u.deg + edge2,\n", - " spherePoint.getRa().asDegrees() * u.deg + edge1/cosd,\n", - " spherePoint.getDec().asDegrees() * u.deg - edge2)\n", + "sq3.polygon = (cutout_ra - Radius/cosd,\n", + " cutout_dec - Radius,\n", + " cutout_ra - Radius/cosd,\n", + " cutout_dec + Radius,\n", + " cutout_ra + Radius/cosd,\n", + " cutout_dec + Radius,\n", + " cutout_ra + Radius/cosd,\n", + " cutout_dec - Radius)\n", "\n", "cutout_bytes = sq3.execute_stream().read()\n", "sq3.raise_if_error()\n", "mem = MemFileManager(len(cutout_bytes))\n", "mem.setData(cutout_bytes, len(cutout_bytes))\n", - "exposure3 = ExposureF(mem)" + "polygon3 = ImageF(mem)" ] }, { @@ -647,17 +720,17 @@ "plt.sca(ax[0])\n", "display1 = afwDisplay.Display(frame=fig)\n", "display1.scale('linear', 'zscale')\n", - "display1.mtv(exposure.image, title='Cutout defined with circle')\n", + "display1.mtv(image, title='Cutout defined with circle')\n", "\n", "plt.sca(ax[1])\n", "display2 = afwDisplay.Display(frame=fig)\n", "display2.scale('linear', 'zscale')\n", - "display2.mtv(exposure2.image, title='Cutout defined with polygon')\n", + "display2.mtv(polygon2, title='Cutout defined with polygon')\n", "\n", "plt.sca(ax[2])\n", "display3 = afwDisplay.Display(frame=fig)\n", "display3.scale('linear', 'zscale')\n", - "display3.mtv(exposure3.image, title='Polygon including cos(dec)')\n", + "display3.mtv(polygon3, title='Polygon including cos(dec)')\n", "plt.tight_layout()\n", "plt.show()" ] @@ -667,7 +740,7 @@ "id": "dfc8a459-e168-40ab-81a5-1148e7536455", "metadata": {}, "source": [ - "> **Figure 3:** A comparison of a cutout generated with the circle function (same as Figure 1; left panel) with a cutout defined using the `polygon` functionality defined using the same edge size (middle panel). The right panel uses `polygon` but accounts for the $cos({\\rm dec})$ term, to replicate the same cutout that was made using circle." + "> **Figure 4:** A comparison of a cutout generated with the circle function (same as Figure 1; left panel) with a cutout defined using the `polygon` functionality defined using the same edge size (middle panel). The right panel uses `polygon` but accounts for the $cos({\\rm dec})$ term, to replicate the same cutout that was made using circle." ] }, { @@ -683,11 +756,20 @@ "id": "99e75a7a-348e-409e-88f6-871b68a4e9dc", "metadata": {}, "source": [ - "## 4. Single extension image cutouts\n", + "## 4. Other image cutout types\n", "\n", - "As mentioned in Section 3, there are three types of cutout services: `cutout-sync-exposure` to return an `ExposureF` type image with all LSST image extensions (demonstrated above); `cutout-sync` which returns just the image extension. The latter will be demonstrated below.\n", + "As mentioned in Section 3, there are three types of cutout services. Demonstrated above is the default `cutout-sync` which returns just the image extension. The cells below demonstrate both the `cutout-sync-exposure` service that returns an `ExposureF` type image with all LSST image extensions (Section 4.1) and the `cutout-sync-mask` service that returns an `MaskedImageF` type containing the LSST image mask (Section 4.2). \n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "a9c8f2cf-f3ca-49af-979e-ee6e8195e034", + "metadata": {}, + "source": [ + "### 4.1 ExposureF cutouts\n", "\n", - "First, reuse the data link `dl_result` that was defined in Section 3 to define a new cutout, this time using `cutout-sync`. Call this soda query `sq4` to differentiate from the earlier ones. Define the cutout as a square subtended by a circle of radius 0.01 degrees, as done in Section 3.\n" + "First, reuse the data link `dl_result` that was defined in Section 3 to define a new cutout, this time using `cutout-sync-exposure` service. Call this soda query `sqE` to differentiate from the earlier ones. Define the cutout as a square subtended by a circle of radius 0.01 degrees, as done in Section 3.\n" ] }, { @@ -697,19 +779,14 @@ "metadata": {}, "outputs": [], "source": [ - "sq4 = SodaQuery.from_resource(dl_result,\n", - " dl_result.get_adhocservice_by_id(\"cutout-sync\"),\n", + "sqE = SodaQuery.from_resource(dl_result,\n", + " dl_result.get_adhocservice_by_id(\"cutout-sync-exposure\"),\n", " session=get_pyvo_auth())\n", "\n", - "spherePoint = geom.SpherePoint(target_ra*geom.degrees, target_dec*geom.degrees)\n", - "Radius = 0.01 * u.deg\n", + "sqE.circle = (cutout_ra, cutout_dec, Radius)\n", "\n", - "sq4.circle = (spherePoint.getRa().asDegrees() * u.deg,\n", - " spherePoint.getDec().asDegrees() * u.deg,\n", - " Radius)\n", - "\n", - "cutout_bytes = sq4.execute_stream().read()\n", - "sq4.raise_if_error()\n" + "cutout_bytes = sqE.execute_stream().read()\n", + "sqE.raise_if_error()\n" ] }, { @@ -733,80 +810,172 @@ }, { "cell_type": "markdown", - "id": "2ae72a88-8002-46b5-8eee-164a2c903314", + "id": "c0eea8e5-ae86-456d-9d67-798557a4a707", "metadata": {}, "source": [ - "Use the `getData` method from the `MemFileManager` class to return and store the bytes associated with the cutout. " + "To enable image display, first store the bytes from `mem` as an LSST `ExposureF` object. The image extension can then be accessed using the `.image` method." ] }, { "cell_type": "code", "execution_count": null, - "id": "ed3f7908-a256-42f5-939f-61f38645c7ca", + "id": "78e1c658-916c-4a73-8c0b-a50d281f4764", "metadata": {}, "outputs": [], "source": [ - "data = mem.getData()" + "exposure = ExposureF(mem)\n" ] }, { "cell_type": "markdown", - "id": "441a49fb-42be-4f37-853e-141419f1da14", + "id": "19815771-fabd-44b1-b359-c1a2f99fcb04", "metadata": {}, "source": [ - "Use the astropy.io fits package to access the bytes and store in an HDUList object. Then, the image extension `image_hdu`, and pixel values `image_data` (a numpy array) can be accessed similar to opening a fits file. " + "The WCS information can be obtained using the getWcs method. " ] }, { "cell_type": "code", "execution_count": null, - "id": "d23a0397-6500-4e36-b956-7917ed6a7173", + "id": "d0d058d6-1cad-4ca4-866b-841baba94220", "metadata": {}, "outputs": [], "source": [ - "hdul = fits.open(io.BytesIO(data))\n", - "hdul.info()\n", + "exposureWCS = WCS(exposure.getWcs().getFitsMetadata())\n", + "print(exposureWCS)" + ] + }, + { + "cell_type": "markdown", + "id": "6a19e7cd-4096-44b9-832f-baa6b6c40482", + "metadata": {}, + "source": [ + "Display the image cutout." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1aa9009a-b493-41ea-bc3c-e0c79adef913", + "metadata": {}, + "outputs": [], + "source": [ + "display = afwDisplay.Display()\n", + "display.scale('asinh', 'zscale')\n", + "display.image(exposure.image)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "ceec7aac-0d64-4397-99df-4450e0098f91", + "metadata": {}, + "source": [ + "> **Figure 5:** The cutout image retrieved in ExposureF format, displayed in pixel coordinates using LSST pipeline tools in grayscale with a scale bar at right.\n", + "\n", + "\n", + "Option to save the entire ExposureF to a fitsfile. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "68c12665-780f-4c47-a16f-e5fd81bc5b77", + "metadata": {}, + "outputs": [], + "source": [ + "# tempdir = os.path.join(os.getenv('HOME'), 'cutouts_temp/')\n", + "# if not os.path.exists(tempdir):\n", + "# os.makedirs(tempdir)\n", + "# print('Created ', tempdir)\n", + "# else:\n", + "# print('Directory already existed: ', tempdir)\n", "\n", - "image_hdu = hdul[1]\n", - "image_data = image_hdu.data\n" + "# sodaCutout = os.path.join(tempdir, 'cutout-circle.fits')\n", + "# with open(sodaCutout, 'bw') as f:\n", + "# f.write(sq.execute_stream().read())" ] }, { "cell_type": "markdown", - "id": "c0eea8e5-ae86-456d-9d67-798557a4a707", + "id": "1d11f7ee-260e-4f52-87fb-e55db25d5d53", "metadata": {}, "source": [ - "To enable image display, first store `image_data` as an LSST `ImageF` object." + "**Note:** `ExposureF` contains multiple Header-Data Units (HDUs) with extensive exposure metadata. To reduce file size, one option is to save only the image, mask, and variance planes, along with the WCS recorded in the Primary HDU. This can be done by writing `exposure.maskedImage`-a three-plane object containing the image, mask, and variance-into a FITS file using `writeFitsImage`." ] }, { "cell_type": "code", "execution_count": null, - "id": "78e1c658-916c-4a73-8c0b-a50d281f4764", + "id": "7b0a245e-3fce-4b35-9ed8-69763b7f1dac", "metadata": {}, "outputs": [], "source": [ - "image = afwImage.ImageF(image_data.astype(\"float32\"))\n" + "# fits.info(sodaCutout)\n", + "\n", + "# sodaCutout_small = os.path.join(tempdir, 'cutout-circle_small.fits')\n", + "# afwDisplay.writeFitsImage(sodaCutout_small,\n", + "# exposure.maskedImage,\n", + "# wcs=exposure.wcs)" ] }, { "cell_type": "markdown", - "id": "6a19e7cd-4096-44b9-832f-baa6b6c40482", + "id": "1fcde65c-8ba1-4ffc-9501-08c01e331e81", "metadata": {}, "source": [ - "Display the image cutout." + "### 4.2 Mask cutouts\n", + "\n", + "Again, reuse the data link `dl_result` that was defined in Section 3 to define a new cutout, this time using `cutout-sync-maskedimage` service. The returned Call this soda query `sqM` to differentiate from the earlier ones. Define the cutout as a square subtended by a circle of radius 0.01 degrees, as done in Section 3.\n" ] }, { "cell_type": "code", "execution_count": null, - "id": "1aa9009a-b493-41ea-bc3c-e0c79adef913", + "id": "0d86fbef-65fe-44d2-bed4-4b6b4ad3a8d9", + "metadata": {}, + "outputs": [], + "source": [ + "sqM = SodaQuery.from_resource(dl_result,\n", + " dl_result.get_adhocservice_by_id(\"cutout-sync-maskedimage\"),\n", + " session=get_pyvo_auth())\n", + "\n", + "sqM.circle = (cutout_ra, cutout_dec, Radius)\n", + "\n", + "cutout_bytes = sqM.execute_stream().read()\n", + "sqM.raise_if_error()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6e9b2366-7051-4766-83a9-2e8e6a6c6b4b", "metadata": {}, "outputs": [], "source": [ + "mem = MemFileManager(len(cutout_bytes))\n", + "mem.setData(cutout_bytes, len(cutout_bytes))\n" + ] + }, + { + "cell_type": "markdown", + "id": "ee886dbe-b264-400d-a7c4-68fa659c739e", + "metadata": {}, + "source": [ + "To enable image display, first store the bytes from `mem` as an LSST `MaskedImageF` object. In addition to the mask, it also contains the image and variance extensions, which can then be accessed using the `.image` or `.variance` attribute." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "178a5016-7490-4410-bda6-8e1885cc03a6", + "metadata": {}, + "outputs": [], + "source": [ + "mask = MaskedImageF(mem)\n", "display = afwDisplay.Display()\n", "display.scale('asinh', 'zscale')\n", - "display.image(image)\n", + "display.image(mask)\n", "plt.show()" ] }, @@ -820,6 +989,8 @@ "id": "acbe9638-e7b7-4229-8301-3f57eb0b0c5b", "metadata": {}, "source": [ + "> **Figure 6:** The cutout image retrieved in MaskedImageF format, displayed in pixel coordinates using LSST pipeline tools in grayscale with a scale bar at right. Mask is plotted on top of the image pixels in colors. \n", + "\n", "## 5. Exercise for the learner\n", "\n", "Reproduce the cutout below, whose center is (ra, dec) = 59.1, -48.8 with 0.06 degrees on a side.\n",