From 57de42c5e9d9cd496dc5660be40e06c6d82f50b5 Mon Sep 17 00:00:00 2001 From: "P. L. Lim" <2090236+pllim@users.noreply.github.com> Date: Fri, 14 Apr 2023 17:36:49 -0400 Subject: [PATCH] EXP: Build test case for advanced aperture photometry. [ci skip] [rtd skip] --- .../concepts/imviz_advanced_aper_phot.ipynb | 403 ++++++++++++++++++ 1 file changed, 403 insertions(+) create mode 100644 notebooks/concepts/imviz_advanced_aper_phot.ipynb diff --git a/notebooks/concepts/imviz_advanced_aper_phot.ipynb b/notebooks/concepts/imviz_advanced_aper_phot.ipynb new file mode 100644 index 0000000000..924edc7575 --- /dev/null +++ b/notebooks/concepts/imviz_advanced_aper_phot.ipynb @@ -0,0 +1,403 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "1f461b4e", + "metadata": {}, + "source": [ + "This concept notebook is an attempt to create a \"worst case scenario\" for aperture photometry in Jdaviz.\n", + "\n", + "See https://github.com/spacetelescope/jdaviz/issues/2139#issuecomment-1507619222" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b946ec7e", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from astropy import coordinates as coord\n", + "from astropy import units as u\n", + "from astropy.modeling import models\n", + "from astropy.nddata import NDData, block_reduce\n", + "from gwcs import coordinate_frames as cf\n", + "from gwcs.wcs import WCS as GWCS\n", + "from photutils.datasets import make_100gaussians_image, make_wcs\n", + "from regions import PixCoord, CirclePixelRegion, RectanglePixelRegion\n", + "from reproject import reproject_interp\n", + "\n", + "from jdaviz import Imviz\n", + "\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6949a5a5-a67e-49a0-a80c-514ac4515098", + "metadata": {}, + "outputs": [], + "source": [ + "imviz = Imviz()\n", + "imviz.show()" + ] + }, + { + "cell_type": "markdown", + "id": "940038df", + "metadata": {}, + "source": [ + "### Image 1: A scene with regular WCS\n", + "\n", + "A hundred Gaussian objects with FITS WCS. This is the data used in `photutils/aperture/tests/test_stats.py`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "431e826e", + "metadata": {}, + "outputs": [], + "source": [ + "data = make_100gaussians_image()\n", + "w = make_wcs(data.shape)\n", + "image_1 = NDData(data, wcs=w)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ea52b644", + "metadata": {}, + "outputs": [], + "source": [ + "imviz.load_data(image_1, data_label=\"fits_wcs\")" + ] + }, + { + "cell_type": "markdown", + "id": "f03325fd-11ac-4e25-a88a-51482c6bfbc4", + "metadata": {}, + "source": [ + "We will use circular apertures from the test case for now so we can check against answers from `photutils`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7cca112b-0cbc-4c66-9948-6846e82cf860", + "metadata": {}, + "outputs": [], + "source": [ + "regions = []\n", + "positions = [(145.1, 168.3), (84.7, 224.1), (48.3, 200.3)]\n", + "for x, y in positions:\n", + " regions.append(CirclePixelRegion(center=PixCoord(x=x, y=y), radius=5))" + ] + }, + { + "cell_type": "markdown", + "id": "ad7b4322-8b51-43b5-ada6-e25f454e0f3c", + "metadata": {}, + "source": [ + "Add one rectangular region too so we can study any assymetrical effects on Subset." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0f18ac72-1424-4820-918f-adf3b29113bf", + "metadata": {}, + "outputs": [], + "source": [ + "regions.append(RectanglePixelRegion(\n", + " center=PixCoord(x=229, y=152), width=17, height=5))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1e1a50f9-66f8-441b-abb7-11091ab51ceb", + "metadata": {}, + "outputs": [], + "source": [ + "imviz.load_regions(regions)" + ] + }, + { + "cell_type": "markdown", + "id": "8e6fc84a-a53d-4de3-ba26-8d6ee4285e02", + "metadata": {}, + "source": [ + "### Image 2: GWCS with rotation and different pixel scale\n", + "\n", + "Now we need an image would result in the Subset being hard to work with. Say, this image would have a different rotation and pixel scale. Offset of the center is also implicitly tested by using a different pixel scale.\n", + "\n", + "*Note: Originally, distortion was also planned but it added a huge amount of computation time to reproject and linking; hence it was abandoned from the test case.*\n", + "\n", + "We start with different pixel scale by down-sampling it by half each dimension." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1185f672-52d4-42c6-b952-7fe587287b07", + "metadata": {}, + "outputs": [], + "source": [ + "# Properties from Image 1\n", + "print(data.shape)\n", + "w.to_header()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b5a18a91-7eb1-483e-ae0a-ab3a78e6596c", + "metadata": {}, + "outputs": [], + "source": [ + "# We want it for the new pixel scale.\n", + "new_shape = tuple((sh // 2) for sh in data.shape)\n", + "new_shape" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6120c794-a00a-45f4-8b49-634e69b7774d", + "metadata": {}, + "outputs": [], + "source": [ + "# https://gwcs.readthedocs.io/en/latest/index.html#getting-started\n", + "\n", + "shift_by_crpix = (models.Shift(-(new_shape[1] - 125 - 0.75) * u.pix) &\n", + " models.Shift(-(new_shape[0] - 75 - 0.75) * u.pix))\n", + "\n", + "matrix = np.array([[-1.3888888888888893e-05, 2.4056261216234408e-05],\n", + " [2.4056261216234408e-05, 1.3888888888888893e-05]])\n", + "rotation = models.AffineTransformation2D(matrix * u.deg,\n", + " translation=[0, 0] * u.deg)\n", + "rotation.input_units_equivalencies = {\"x\": u.pixel_scale(2 * (u.deg/u.pix)),\n", + " \"y\": u.pixel_scale(2 * (u.deg/u.pix))}\n", + "rotation.inverse = models.AffineTransformation2D(np.linalg.inv(matrix) * u.pix,\n", + " translation=[0, 0] * u.pix)\n", + "rotation.inverse.input_units_equivalencies = {\"x\": u.pixel_scale(0.5 * (u.pix / u.deg)),\n", + " \"y\": u.pixel_scale(0.5 * (u.pix / u.deg))}\n", + "\n", + "tan = models.Pix2Sky_TAN()\n", + "celestial_rotation = models.RotateNative2Celestial(\n", + " 197.8925 * u.deg, -1.36555556 * u.deg, 180 * u.deg)\n", + "\n", + "det2sky = shift_by_crpix | rotation | tan | celestial_rotation\n", + "det2sky.name = \"linear_transform\"\n", + "\n", + "detector_frame = cf.Frame2D(name=\"detector\", axes_names=(\"x\", \"y\"),\n", + " unit=(u.pix, u.pix))\n", + "sky_frame = cf.CelestialFrame(reference_frame=coord.ICRS(), name='icrs',\n", + " unit=(u.deg, u.deg))\n", + "\n", + "pipeline = [(detector_frame, det2sky),\n", + " (sky_frame, None)]\n", + "gw = GWCS(pipeline)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bac3b7c2-a9e0-4b30-a092-826e00078862", + "metadata": {}, + "outputs": [], + "source": [ + "data_sm = block_reduce(data, 2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "874e9cef-85ea-4055-b1d5-3802efcba7a8", + "metadata": {}, + "outputs": [], + "source": [ + "image_2 = NDData(data_sm, wcs=gw)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "05832fda-22b4-42de-98fb-254ef396b03d", + "metadata": {}, + "outputs": [], + "source": [ + "imviz.load_data(image_2, data_label=\"gwcs\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8c2793e8-65d4-46ed-a946-347767d1309d", + "metadata": {}, + "outputs": [], + "source": [ + "# imviz.link_data(link_type='wcs')" + ] + }, + { + "cell_type": "markdown", + "id": "987c96f0-ba8b-42b1-9ffd-6f50db3f190c", + "metadata": {}, + "source": [ + "### Image 2.1: Add rotation to downsampled image" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "91dc1847-1ae1-4a51-aa44-feac12f4c290", + "metadata": {}, + "outputs": [], + "source": [ + "# Like above but with extra rotation\n", + "celestial_rotation = models.RotateNative2Celestial(\n", + " 197.8925 * u.deg, -1.36555556 * u.deg, 190 * u.deg)\n", + "\n", + "det2sky = shift_by_crpix | rotation | tan | celestial_rotation\n", + "det2sky.name = \"linear_transform\"\n", + "\n", + "pipeline = [(detector_frame, det2sky),\n", + " (sky_frame, None)]\n", + "gw_rot = GWCS(pipeline)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "43797cca-c6e1-4bc0-ad66-46ba5b16bcfe", + "metadata": {}, + "outputs": [], + "source": [ + "# This takes much longer when GWCS has distortion.\n", + "data_rot, _ = reproject_interp(image_2, gw_rot, shape_out=data_sm.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1b52e436-2cd3-4a42-b7a8-d489a6fb06ad", + "metadata": {}, + "outputs": [], + "source": [ + "image_2_1 = NDData(data_rot, wcs=gw_rot)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "42cadf07-99cc-4816-b6c3-43847c73478e", + "metadata": {}, + "outputs": [], + "source": [ + "imviz.load_data(image_2_1, data_label=\"gwcs_rot\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8d78fcfe-d224-410a-a9d8-7ff84a3bd660", + "metadata": {}, + "outputs": [], + "source": [ + "imviz.link_data(link_type='wcs')" + ] + }, + { + "cell_type": "markdown", + "id": "fe16d7c9-5b41-45fe-b7a5-3f6c1e91ec95", + "metadata": {}, + "source": [ + "## How is glue projecting the rectangle?\n", + "\n", + "How is `glue` really projecting the rectangle for all the data?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "69549536-96f5-43ea-afc4-aaa6f79723ba", + "metadata": {}, + "outputs": [], + "source": [ + "rect_grp = imviz.app.data_collection.subset_groups[3]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5958db68-2d75-49da-bac5-4a1eefb59ded", + "metadata": {}, + "outputs": [], + "source": [ + "fig, axs = plt.subplots(1, 3)\n", + "axs[0].imshow(rect_grp.subsets[0].to_mask(), origin='lower')\n", + "axs[0].set_xlim(210, 250)\n", + "axs[0].set_ylim(140, 160)\n", + "axs[1].imshow(rect_grp.subsets[1].to_mask(), origin='lower')\n", + "axs[1].set_xlim(105, 125)\n", + "axs[1].set_ylim(70, 80)\n", + "axs[2].imshow(rect_grp.subsets[2].to_mask(), origin='lower')\n", + "axs[2].set_xlim(105, 125)\n", + "axs[2].set_ylim(70, 80);" + ] + }, + { + "cell_type": "markdown", + "id": "2071b1be-ee74-41c3-bb36-02eb1eb4f0a8", + "metadata": {}, + "source": [ + "Even when unrotated, a different pixel scale does affect the projected mask dimension. When rotated, it is not even a rectangle anymore; furthermore, even though this projection follows WCS linking, the array from `to_mask()` does not account for rotation that we see in the viewer." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5efe10a8-4e83-4c88-92ea-cc7ad4b42538", + "metadata": {}, + "outputs": [], + "source": [ + "# TODO: Now how do we use this test case?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "192e6319-8af3-4450-980f-61fca40f711f", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "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.10.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}