{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Tutorial / How to use\n", "\n", "In this tutorial we create a (simplified) synthetic galaxy image from scratch, along with its associated segmentation map, and then run the statmorph code on it.\n", "\n", "If you already have a real astronomical image and segmentation map to work with, jump to [Running statmorph](#Running-statmorph).\n", "\n", "\n", "### Setting up\n", "\n", "We import some Python packages first. If you are missing any of these, please see the the installation instructions." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from astropy.visualization import simple_norm\n", "from astropy.modeling.models import Sersic2D\n", "from astropy.convolution import convolve, Gaussian2DKernel\n", "from photutils.segmentation import detect_threshold, detect_sources\n", "import time\n", "import statmorph\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Creating a model galaxy image**\n", "\n", "We assume that the image size is 240x240 pixels and that the \"true\" light distribution is described by a 2D Sersic model with the following parameters:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ny, nx = 240, 240\n", "y, x = np.mgrid[0:ny, 0:nx]\n", "sersic_model = Sersic2D(\n", " amplitude=1, r_eff=20, n=2.5, x_0=120.5, y_0=96.5,\n", " ellip=0.5, theta=0.5)\n", "image = sersic_model(x, y)\n", "plt.imshow(image, cmap='gray', origin='lower',\n", " norm=simple_norm(image, stretch='log', log_a=10000))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Convolving with a PSF**\n", "\n", "In practice, every astronomical image is the convolution of a \"true\" image with a point spread function (PSF), which depends on the optics of the telescope, atmospheric conditions, etc. Here we assume that the PSF is a simple 2D Gaussian kernel with a standard deviation of 2 pixels:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "kernel = Gaussian2DKernel(2)\n", "kernel.normalize() # make sure kernel adds up to 1\n", "psf = kernel.array # we only need the numpy array\n", "plt.imshow(psf, origin='lower', cmap='gray')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we convolve the image with the PSF." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "image = convolve(image, psf)\n", "plt.imshow(image, cmap='gray', origin='lower',\n", " norm=simple_norm(image, stretch='log', log_a=10000))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Applying shot noise**\n", "\n", "One source of noise in astronomical images originates from the Poisson statistics of the number of electrons recorded by each pixel. We can model this effect by introducing a *gain* parameter, a scalar that can be multiplied by the science image to obtain the number of electrons per pixel.\n", "\n", "For the sake of this example, we choose a very large gain value, so that shot noise becomes almost negligible (10^5 electrons/pixel at the effective radius, where we had defined an amplitude of 1.0 in arbitrary units). The resulting image after applying shot noise looks very similar to the one from the previous step and is not shown." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "np.random.seed(3)\n", "gain = 1e5\n", "image = np.random.poisson(image * gain) / gain" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Applying background noise**\n", "\n", "Apart from shot noise, astronomical images have a sky background noise component, which we here model with a uniform Gaussian distribution centered at zero (since the image is background-subtracted).\n", "\n", "We assume, somewhat optimistically, that the signal-to-noise ratio (S/N) per pixel is 100 at the effective radius (where we had defined the Sersic model amplitude as 1.0)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "snp = 100.0\n", "sky_sigma = 1.0 / snp\n", "image += sky_sigma * np.random.standard_normal(size=(ny, nx))\n", "plt.imshow(image, cmap='gray', origin='lower',\n", " norm=simple_norm(image, stretch='log', log_a=10000))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Gain and weight maps**\n", "\n", "Note that statmorph will ask for one of two input arguments: (1) a weight map, which is a 2D array (of the same size as the input image) representing one standard deviation at each pixel value, or (2) the gain, which was described above. In the latter case, the gain parameter is used internally by statmorph (along with an automatic estimation of the sky background noise) to calculate the weight map." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Creating a segmentation map**\n", "\n", "Besides the image itself and the weight map/gain, the only other *required* argument is the segmentation map, which labels the pixels belonging to different sources. It is usually generated with specialized tools such as SExtractor or photutils. Here we use the latter to create a simplified segmentation map, where we detect sources that lie above a 1.5-sigma detection threshold.\n", "\n", "Note that the detection stage (but, importantly, not the threshold calculation) is carried out on a \"convolved\" version of the image. This is done in order to smooth out small-scale noise and thus ensure that the shapes of the detected segments are reasonably smooth." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "threshold = detect_threshold(image, 1.5)\n", "npixels = 5 # minimum number of connected pixels\n", "convolved_image = convolve(image, psf)\n", "segmap = detect_sources(convolved_image, threshold, npixels)\n", "plt.imshow(segmap, origin='lower', cmap='gray')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this particular case, the obtained segmap has only two values: 0 for the background (as should always be the case) and 1 for the only labeled source. However, statmorph is designed to process all the sources labeled by a segmentation map, which makes it applicable to large mosaic images." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Running statmorph\n", "\n", "Now that we have all the required data, we are ready to measure the morphology of the source just created. Note that we include the PSF as a keyword argument, which results in more correct Sersic profile fits.\n", "\n", "Also note that we do not attempt to fit a *double* Sersic model, which would be degenerate in this particular case (the two components would be identical and their relative amplitudes would be unconstrained). For a demonstration of statmorph's double Sersic fitting functionality, see the [Double 2D Sersic example](./doublesersic.html)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "start = time.time()\n", "source_morphs = statmorph.source_morphology(\n", " image, segmap, gain=gain, psf=psf)\n", "print('Time: %g s.' % (time.time() - start))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Examining the output\n", "\n", "In general, `source_morphs` is a list of objects corresponding to each labeled source in the image. Here we focus on the first (and only) labeled source:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "morph = source_morphs[0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we print some of the morphological properties just calculated:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print('BASIC MEASUREMENTS (NON-PARAMETRIC)')\n", "print('xc_centroid =', morph.xc_centroid)\n", "print('yc_centroid =', morph.yc_centroid)\n", "print('ellipticity_centroid =', morph.ellipticity_centroid)\n", "print('elongation_centroid =', morph.elongation_centroid)\n", "print('orientation_centroid =', morph.orientation_centroid)\n", "print('xc_asymmetry =', morph.xc_asymmetry)\n", "print('yc_asymmetry =', morph.yc_asymmetry)\n", "print('ellipticity_asymmetry =', morph.ellipticity_asymmetry)\n", "print('elongation_asymmetry =', morph.elongation_asymmetry)\n", "print('orientation_asymmetry =', morph.orientation_asymmetry)\n", "print('rpetro_circ =', morph.rpetro_circ)\n", "print('rpetro_ellip =', morph.rpetro_ellip)\n", "print('rhalf_circ =', morph.rhalf_circ)\n", "print('rhalf_ellip =', morph.rhalf_ellip)\n", "print('r20 =', morph.r20)\n", "print('r80 =', morph.r80)\n", "print('Gini =', morph.gini)\n", "print('M20 =', morph.m20)\n", "print('F(G, M20) =', morph.gini_m20_bulge)\n", "print('S(G, M20) =', morph.gini_m20_merger)\n", "print('sn_per_pixel =', morph.sn_per_pixel)\n", "print('C =', morph.concentration)\n", "print('A =', morph.asymmetry)\n", "print('S =', morph.smoothness)\n", "print()\n", "print('SERSIC MODEL')\n", "print('sersic_amplitude =', morph.sersic_amplitude)\n", "print('sersic_rhalf =', morph.sersic_rhalf)\n", "print('sersic_n =', morph.sersic_n)\n", "print('sersic_xc =', morph.sersic_xc)\n", "print('sersic_yc =', morph.sersic_yc)\n", "print('sersic_ellip =', morph.sersic_ellip)\n", "print('sersic_theta =', morph.sersic_theta)\n", "print('sersic_chi2_dof =', morph.sersic_chi2_dof)\n", "print()\n", "print('OTHER')\n", "print('sky_mean =', morph.sky_mean)\n", "print('sky_median =', morph.sky_median)\n", "print('sky_sigma =', morph.sky_sigma)\n", "print('flag =', morph.flag)\n", "print('flag_sersic =', morph.flag_sersic)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the fitted Sersic model is in very good agreement with the \"true\" Sersic model that we originally defined (n = 2.5, rhalf = 20, etc.) and that the reduced chi-squared statistic (sersic_chi2_dof) is close to 1, indicating a good fit without overfitting. However, such good agreement tends to deteriorate somewhat at higher noise levels, and one has to keep in mind that not all galaxies are well described by Sersic profiles.\n", "\n", "Other morphological measurements that are more general and robust to noise, which are also calculated by statmorph, include the Gini-M20 (Lotz et al. 2004), CAS (Conselice 2003) and MID (Freeman et al. 2013) statistics, as well as the outer asymmetry (Wen et al. 2014) and shape asymmetry (Pawlik et al. 2016).\n", "\n", "Also note that statmorph returns two quality flags:\n", "\n", "1. ``flag`` : indicates the quality of the basic morphological measurements, taking one of the following values: 0 (good), 1 (suspect), 2 (bad), or 4 (catastrophic). More details can be found [here](../description.html#output).\n", "\n", "2. ``flag_sersic`` : indicates the quality of the Sersic fit, also taking values of 0 (good), 1 (suspect), 2 (bad), or 4 (catastrophic).\n", "\n", "In general, ``flag <= 1`` should always be enforced, while ``flag_sersic <= 1`` should only be used when one is interested in the Sersic fits (which might fail for merging galaxies and other \"irregular\" objects)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Visualizing the morphological measurements**\n", "\n", "For convenience, statmorph includes a ``make_figure`` function that can be used to visualize some of the morphological measurements. This creates a multi-panel figure analogous to Fig. 4 from Rodriguez-Gomez et al. (2019)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from statmorph.utils.image_diagnostics import make_figure\n", "fig = make_figure(morph)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig.savefig('tutorial.png', dpi=150)\n", "plt.close(fig)" ] } ], "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.11.6" } }, "nbformat": 4, "nbformat_minor": 4 }