Merge pull request #52 from tddesjardins/galshapes

Galaxy shapes update

content/notebooks/measuring_galaxy_shapes/measuring_galaxy_shapes.ipynb

@@ -41,6 +41,7 @@
  "- *astropy.table* for manipulating tabular data\n",
  "- *astropy.units* for working with units\n",
  "- *astropy.visualization* for normalizing images for plotting\n",
+ "- *coord* for working with angles\n",
  "- *matplotlib.pyplot* for plotting and visualizing data\n",
  "- *galsim* to measure source momens and generate Sérsic profiles\n",
  "- *numpy* for array manipulation and mathematical operations\n",
@@ -67,6 +68,7 @@
  "from astropy.table import Table\n",
  "import astropy.units as u\n",
  "from astropy.visualization import simple_norm\n",
+ "import coord\n",
  "import matplotlib.pyplot as plt\n",
  "import galsim\n",
  "from galsim.roman import collecting_area\n",
@@ -92,13 +94,17 @@
  "## Introduction\n",
  "The main goal of this notebook is to illustrate a typical use case of Roman images, which is performing shape measurements of astronomical sources. \n",
  "\n",
- "We are going to perform two sets of measurements. First, we rely on [`galsim`](https://galsim-developers.github.io/) to perform ellipticity measurements. In particular, we use the\n",
+ "We are going to perform two sets of measurements. First, we rely on [GalSim](https://galsim-developers.github.io/) to perform ellipticity measurements. In particular, we use the\n",
  "[REGAUSS method in its HSM module](https://galsim-developers.github.io/GalSim/_build/html/hsm.html) ([Hirata & Seljak 2003](https://ui.adsabs.harvard.edu/abs/2003MNRAS.343..459H/abstract), [Mandelbaum et al. 2005](https://ui.adsabs.harvard.edu/abs/2005MNRAS.361.1287M/abstract)). Second, we fit a Sérsic model to a galaxy cutout.\n",
  "\n",
- "This notebook relies on ASDF file manipulations. For additional information about these files\n",
- "and how to use them, please check our Working with ASDF notebook tutorial.\n",
+ "In this notebook, we will be building on previous tutorials. We recommend consulting the following tutorials prior to this one:\n",
+ "- Data Access and Discovery\n",
+ "- Working with ASDF\n",
+ "- Data Visualization\n",
+ "- WebbPSF\n",
+ "- Synphot\n",
  "\n",
- "### Defining terms\n",
+ "### Defining Terms\n",
  "\n",
  "- ASDF: Advanced Scientific Data Format — the data format for the Roman Wide Field Instrument"
  ]
@@ -118,7 +124,7 @@
  }
  },
  "source": [
- "## Loading data\n",
+ "## Loading Data\n",
  "The first step of the analysis is to read the Roman WFI image data, which are stored in ASDF format. For this example, we start with a calibrated Level 2 (L2) simulated image created with [Roman I-Sim](https://romanisim.readthedocs.io). For more information about Roman's data products check the [Roman User Documentation](roman-docs.stsci.edu). We use `roman_datamodels` to open the simulated image, which we stream into memory from an S3 bucket:"
  ]
  },
@@ -207,10 +213,10 @@
  "cell_type": "markdown",
  "metadata": {},
  "source": [
- "## Creating cutouts and retrieving WebbPSF's PSF\n",
+ "## Creating Cutouts and Retrieving the PSF\n",
  "\n",
- "The positions are in sky coordinates. In order to make the cutouts we need the pixel coordinates in the image. In order to do that, we use the input World Coordinate System (WCS).\n",
- "In the case of Roman data, these correspond to GWCS objects. For more information about GWCS please check the documentation [here](https://gwcs.readthedocs.io/)."
+ "The positions are in sky coordinates. In order to make the cutouts we need the pixel coordinates in the image. To do that transformation, we will use the generalized World Coordinate System (GWCS) object in the simulated file's metadata.\n",
+ "For more information about GWCS please check the documentation [here](https://gwcs.readthedocs.io/)."
  ]
  },
  {
@@ -267,7 +273,7 @@
  "cell_type": "markdown",
  "metadata": {},
  "source": [
- "Now let's make our subselection of our catalog:"
+ "Now let's trim our catalog based on our magnitude cuts and selection of extended sources:"
  ]
  },
  {
@@ -352,7 +358,7 @@
  "cell_type": "markdown",
  "metadata": {},
  "source": [
- "Now, before the next step of creating a cutout, we need to make an estimate of the sky background level. Calibrated WFI images are not background subtracted by the pipeline, so we need to do this ourselves. There are various ways to accomplish this, but for our example we'll use the sigma_clipped median of the image:"
+ "Now, before the next step of creating a cutout, we need to make an estimate of the sky background level. Calibrated WFI images are not background subtracted by the pipeline, so we need to do this ourselves. There are various ways to accomplish this, but for our example we'll use the sigma_clipped median of the image. This background estimate is simplistic, but effective for our demonstration."
  ]
  },
  {
@@ -434,11 +440,11 @@
  "cell_type": "markdown",
  "metadata": {},
  "source": [
- "As we can see, the background estimate seems to be reasonable. We likely could improve this by doing a local background estimate rather than using the whole image. A significant change in our background subtraction would affect the results of our later model fitting, and so background estimation should be treated with care.\n",
+ "As we can see, the background estimate seems to be reasonable. We likely could improve this by doing a local background estimate rather than using the whole image. In real, on-orbit observations, a 2-D local background model would likely provide the highest fidelity background estimate. A significant change in our background subtraction would affect the results of our later model fitting, and so background estimation should be treated with care.\n",
  "\n",
- "## Estimating the source's moments\n",
+ "## Estimating the Source Moments\n",
  "\n",
- "Now that we have our cutout saved as a `galsim.Image` object we can just use the `HSM` module to estimate the moments:"
+ "Now that we have our cutout saved as a `galsim.Image` object we can just use the `hsm` module to estimate the moments:"
  ]
  },
  {
@@ -455,7 +461,30 @@
  "cell_type": "markdown",
  "metadata": {},
  "source": [
- "The `shape` variable contains the source moments estimated from the best-fit elliptical Gaussian. Since we have the PSF model we can also try to estimate the PSF-corrected galaxy shear:"
+ "The `shape` variable contains the source moments estimated from the best-fit elliptical Gaussian. For example, we can see in the printed result that the `.observed_shape` attribute contains a `galsim.Shear` object. From that object, we can retrieve shape parameters such as the position angle (beta) and the minor-to-major axis ratio (q):"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Get the position angle (beta) and wrap it to the range [0, 2pi) radians and convert to degrees\n",
+ "beta = shape.observed_shape.beta.wrap(center = 180 * coord.degrees).deg\n",
+ "\n",
+ "# Get the minor-to-major axis ratio (q)\n",
+ "q = shape.observed_shape.q\n",
+ "\n",
+ "print(f'position angle: {beta:.5f} deg')\n",
+ "print(f'minor-to-major axis ratio: {q:.5f}')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Since we have the PSF model we can also try to estimate the PSF-corrected galaxy shear:"
  ]
  },
  {
@@ -472,18 +501,43 @@
  "cell_type": "markdown",
  "metadata": {},
  "source": [
- "## Fitting a source to a Sérsic profile"
+ "We can once again extract information about the position angle and minor-to-major axis ratio as we did before:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "beta = shape2.observed_shape.beta.wrap(center = 180 * coord.degrees).deg\n",
+ "q = shape2.observed_shape.q\n",
+ "print(f'position angle: {beta:.5f} deg')\n",
+ "print(f'minor-to-major axis ratio: {q:.5f}')"
  ]
  },
  {
  "cell_type": "markdown",
  "metadata": {},
  "source": [
- "Another typical morphological analysis consists on fitting the source with a [Sérsic profile](https://en.wikipedia.org/wiki/S%C3%A9rsic_profile) ([Sérsic 1963](https://ui.adsabs.harvard.edu/abs/1963BAAA....6...41S/abstract)). We are going to rely on the Sérsic implementation in the `galsim` package instead of doing our own.\n",
+ "In this case, we see that the values are the same with and without the PSF-correction.\n",
+ "\n",
+ "## Fitting a Sérsic Profile to a Source"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Another typical morphological analysis consists of fitting a source with a [Sérsic profile](https://en.wikipedia.org/wiki/S%C3%A9rsic_profile) ([Sérsic 1963](https://ui.adsabs.harvard.edu/abs/1963BAAA....6...41S/abstract)). We are going to rely on the Sérsic implementation in the `galsim` package instead of doing our own.\n",
  "\n",
  "In our example we will ignore the pixel mask (data quality array), but will use the error image.\n",
  "\n",
- "First, we need to convert our source flux units to photons / s / cm$^2$. The catalog fluxes are stored as maggies, which are linear flux units such that $m_\\mathrm{AB} = -2.5\\times\\log_{10}(f)$ where $m_\\mathrm{AB}$ is the AB magnitude and $f$ is the flux in maggies. As Roman I-Sim, which was used to generate the catalog and simulate the image, does not create sources with realistic colors, let's make a simplifying assumption that our source has a spectral energy distribution (SED) that is flat in frequency space, at least over the bandpass of interest (F106). We will use synphot to synthetically observe the source spectrum, normalized to the AB magnitude from the catalog, and integrate over the bandpass:"
+ "First, we need to convert our source flux units to photons / s / cm$^2$. The catalog fluxes are stored as maggies, which are linear flux units such that:\n",
+ "\n",
+ "$m_\\mathrm{AB} = -2.5\\times\\log_{10}(f),$ \n",
+ "\n",
+ "where $m_\\mathrm{AB}$ is the AB magnitude and $f$ is the flux in maggies. As Roman I-Sim, which was used to generate the catalog and simulate the image, does not create sources with realistic colors, let's make a simplifying assumption that our source has a spectral energy distribution (SED) that is flat in frequency space, at least over the bandpass of interest (F106). We will use the `synphot` package to synthetically observe the source spectrum, normalized to the AB magnitude from the catalog, and integrate over the bandpass:"
  ]
  },
  {
@@ -508,22 +562,30 @@
  "# Angstroms. Thus, the integral of the Observation over the wavelengths will\n",
  "# give units of photons / s / cm^2.\n",
  "waves = obs.binset\n",
- "flux = np.trapz(obs(waves), x=waves).value"
+ "flux = np.trapz(obs(waves), x=waves).value\n",
+ "\n",
+ "print(f'Source flux: {flux:.5f} photon / s / cm^2')"
  ]
  },
  {
  "cell_type": "markdown",
  "metadata": {},
  "source": [
- "Then we create our model as a callable function for fitting later. The variables are:\n",
+ "Note that we could have instead determined a scaling factor to apply to any input flux from the catalog by normalizing the spectrum to an amplitude of 1 maggie. Synphot does not allow us to use the `astropy.units.maggy` unit, but we can just as effectively use the definition of AB magnitudes to write this in units of Jankys:\n",
+ "\n",
+ "$m_\\mathrm{AB} = -2.5 \\times log_{10}(f_\\mathrm{Jy}) + 8.90,$\n",
+ "\n",
+ "and thus when $m_\\mathrm{AB} = 0$ mag, which is defined when $f = 1$ maggie, we find that $f_\\mathrm{Jy} = 10^{(8.90 / 2.5)} \\approx 3630.78055~\\mathrm{Jy}$. Determining such a scaling factor would be useful if we wanted to perform our shape analysis at scale.\n",
+ "\n",
+ "Then we create our model as a callable function for later fitting. The variables are:\n",
  "* n = Sérsic index\n",
  "* hlr = half-light radius in arcseconds\n",
  "* pa = position angle in degrees\n",
  "* x0 = shift in x pixels to center of the source\n",
  "* y0 = shift in y pixels to center of the source\n",
  "* q = minor-to-major axis ratio\n",
  "\n",
- "We will set the gain and area options in the `drawImage()` method on the `galsim.Sersic` object we create. This ensures that the model of our galaxy can return pixels with values like our observed image that has units of DN / s."
+ "We will set the gain and area options in the `drawImage()` method on the `galsim.Sersic` object we create. This ensures that the model of our galaxy can return pixels with values like our observed image that has units of DN / s. We set the gain to 1.8, which is a reasonable value for a general WFI detector, and we use the obscured collecting area from the `galsim.roman` module."
  ]
  },
  {
@@ -539,7 +601,7 @@
  " offset = galsim.PositionD(x=x0, y=y0) # potentially shift a bit the profile\n",
  " ser = galsim.convolve.Convolve(ser, psf_obj) # add PSF\n",
  " img = galsim.ImageD(nx, ny, scale=0.11) \n",
- " ser.drawImage(image=img, offset=offset, scale=0.11, gain=2, area=np.pi * (237 / 2)**2)\n",
+ " ser.drawImage(image=img, offset=offset, scale=0.11, gain=1.8, area=galsim.roman.collecting_area)\n",
  " \n",
  " return img.array.flatten()"
  ]
@@ -568,7 +630,7 @@
  "cell_type": "markdown",
  "metadata": {},
  "source": [
- "Notice that we used the syntax `*p0` above. This is the same as explicitly writing out `p0[0], p0[1], ..., p0[6]` and is a convenient way in Python to unpack all of the values in a list in order as inputs to function's positional arguments.\n",
+ "Notice that we used the syntax `*p0` above. This is the same as explicitly writing out `p0[0], p0[1], ..., p0[6]` and is a convenient way in Python to unpack all of the values in a list in order as inputs to a function's positional arguments.\n",
  "\n",
  "We also obtain the error array in the location of the cutout:"
  ]
@@ -626,6 +688,8 @@
  "cell_type": "markdown",
  "metadata": {},
  "source": [
+ "Comparing the input and best-fit parameters, we see good agreement.\n",
+ "\n",
  "Let's create an image of our best-fit model and visualize the original cutout, initial input model, best-fit output model, and fit residuals:"
  ]
  },
@@ -726,7 +790,41 @@
  "cell_type": "markdown",
  "metadata": {},
  "source": [
- "As we can see, the best-fitting model is quite similar to the previous one with the exceptions of the minor-to-major axis ratio and the position angle. This means that these parameters do not greatly influence the model and/or are not well-constrained by the data of our example. In fact, on closer examination, we can see that the position angle and minor-to-major axis ratio in our earlier fit is identical to the input values we obtained from the catalog. This demonstrates the difficulty of measuring these values in observational data.\n",
+ "As we can see, the best-fitting model is quite similar to the previous one with the exceptions of the flux, minor-to-major axis ratio, and the position angle. In fact, on closer examination, we can see that the position angle and minor-to-major axis ratio in our earlier fit is identical to the input values we obtained from the catalog.\n",
+ "\n",
+ "Let's try again, but this time let's use the shear estimate shape measurements we previously performed (recall that beta and q are the position angle and minor-to-major axis ratio, respectively, from the PSF-corrected shear estimate in the `shape2` variable above): "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "p1 = [1., 0.44, flux, beta, 0., 0., q]\n",
+ "\n",
+ "# Create the model with our truth parameters\n",
+ "fid_model_new = sersic_mod(cutout, *p1)\n",
+ "\n",
+ "pout_new, pcov_new = curve_fit(sersic_mod, cutout, cutout.array.flatten(), sigma=err_img,\n",
+ " p0=p1, bounds=([0.3, 0.01, 1e-11, 0, -size/2, -size/2, 0], [6.2, 4.0, 0.1, 360, size/2, size/2, 1]))\n",
+ "\n",
+ "labels = ('n', 'hlr', 'flux', 'pa', 'x0', 'y0', 'q')\n",
+ "\n",
+ "for i, _ in enumerate(p0):\n",
+ " print(f'{labels[i]}:')\n",
+ " print(f'\\tCatalog (p0): {p0[i]:.5f}')\n",
+ " print(f'\\tNew inputs (p1): {p1[i]:.5f}')\n",
+ " print('')\n",
+ " print(f'\\tOld fit: {pout[i]:.5f}')\n",
+ " print(f'\\tNew fit: {pout_new[i]:.5f}')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Much better! Now the best-fit flux agrees well with the catalog value. We also see that we got similar values of the other best-fit parameters other than the position angle and minor-to-major axis ratio. It seems that these two parameters are not well constrained by the data using the Sérsic model fit.\n",
  "\n",
  "For completeness, let's also re-plot the cutout, models, and fit residuals with out new inputs and best fit:"
  ]
@@ -794,7 +892,7 @@
  "cell_type": "markdown",
  "metadata": {},
  "source": [
- "Comparing these values with our earlier fit, we find that the median and standard deviation of the residuals of both fits are identical to within four decimal places despite differences in the best-fit model parameters. "
+ "Comparing these values with our earlier fit, we find that the median and standard deviation of the residuals of both fits are identical to within four decimal places despite differences in the input and best-fit model parameters. "
  ]
  },
  {
@@ -805,8 +903,9 @@
  "\n",
  "- [Roman User Documentation (RDox)](https://roman-docs.stsci.edu/)\n",
  "- [Roman Notebooks](https://github.com/spacetelescope/roman_datamodels)\n",
- "- [`romanisim` documentation](https://romanisim.readthedocs.io/)\n",
- "- [`webbpsf` documentation](https://webbpsf.readthedocs.io/)"
+ "- [Roman I-Sim documentation](https://romanisim.readthedocs.io/)\n",
+ "- [WebbPSF documentation](https://webbpsf.readthedocs.io/)\n",
+ "- [GalSim HSM shape fitting documentation](https://galsim-developers.github.io/GalSim/_build/html/hsm.html)"
  ]
  },
  {
@@ -817,7 +916,7 @@
  }
  },
  "source": [
- "## About this notebook\n",
+ "## About This Notebook\n",
  "\n",
  "**Author:** Javier Sánchez, Amethyst Barnes, Ami Choi, Tyler Desjardins \n",
  "**Updated On:** 2024-09-24"