diff --git a/.gitignore b/.gitignore index e493f4b..2ee9f9a 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ .pyc .ipynb_checkpoints/ *.nbconvert.ipynb -photutils/dev/* \ No newline at end of file +photutils/dev/ + diff --git a/photutils/02_photutils_source_detection.ipynb b/photutils/02_photutils_source_detection.ipynb new file mode 100644 index 0000000..ead2021 --- /dev/null +++ b/photutils/02_photutils_source_detection.ipynb @@ -0,0 +1,1100 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "\n", + "[\n", + "\n", + "](http://photutils.readthedocs.io/en/stable/index.html)\n", + "\n", + "\n", + "# Source Detection with `photutils`\n", + "---" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##### What is source detection?\n", + "In order to do photometry on astronomical image data, one must first determine the locations of the sources in the image. Source detection methods find sources algorithmically, by looking for regions of an image where the signal from a source is statistically higher than the signal from background noise. Some algorithms search for sources whose profiles match specific data models, such as a 2-D Gaussian, while others simply look for local maxima.\n", + "\n", + "The `photutils` package provides a variety of tools that use different detection algorithms to locate sources in an image.\n", + "\n", + "##### What does this tutorial include?\n", + "This tutorial covers different tools for source detection with `photutils`, including the following methods:\n", + "* Source Detection with `DAOStarFinder`\n", + "* Source Detection with `IRAFStarFinder`\n", + "* Source Detection with `find_peaks`\n", + "* Image Segmentation\n", + "\n", + "##### Which data are used in this tutorial?\n", + "We will be manipulating Hubble eXtreme Deep Field (XDF) data, which was collected using the Advanced Camera for Surveys (ACS) on Hubble between 2002 and 2012. The image we use here is the result of 1.8 million seconds (500 hours!) of exposure time, and includes some of the faintest and most distant galaxies that have ever been observed. \n", + "\n", + "##### The methods demonstrated here are available in narrative form within the `photutils.detection` [documentation](http://photutils.readthedocs.io/en/stable/detection.html) and `photutils.segmentation` [documentation](http://photutils.readthedocs.io/en/stable/segmentation.html)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "\n", + "**Important:** Before proceeding, please be sure to update your versions of `astropy`, `matplotlib`, and `photutils`, or this notebook may not work properly. Or, if you don't want to handle packages individually, you can always use (and keep updated!) the [AstroConda](https://astroconda.readthedocs.io) distribution.\n", + "\n", + "
\n", + "\n", + "---" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Import necessary packages" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, let's import packages that we will use to perform arithmetic functions and visualize data:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from astropy.io import fits\n", + "import astropy.units as u\n", + "from astropy.nddata import CCDData\n", + "from astropy.stats import sigma_clipped_stats, SigmaClip\n", + "from astropy.visualization import ImageNormalize, LogStretch\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib.ticker import LogLocator\n", + "import numpy as np\n", + "from photutils.background import Background2D, MeanBackground\n", + "\n", + "# Show plots in the notebook\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's also define some `matplotlib` parameters, such as title font size and the dpi, to make sure our plots look nice. To make it quick, we'll do this by loading a [style file shared with the other photutils tutorials](photutils_notebook_style.mplstyle) into `pyplot`. We will use this style file for all the notebook tutorials. (See [here](https://matplotlib.org/users/customizing.html) to learn more about customizing `matplotlib`.)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.style.use('photutils_notebook_style.mplstyle')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Retrieve data" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As described in the introduction, we will be using Hubble eXtreme Deep Field (XDF) data. Since this file is too large to store on GitHub, we will just use `astropy` to directly download the file from the STScI archive: https://archive.stsci.edu/prepds/xdf/ \n", + "\n", + "(Generally, the best package for web queries of astronomical data is [Astroquery](https://astroquery.readthedocs.io/en/latest/); however, the dataset we are using is a High Level Science Product (HLSP) and thus is not located within a catalog that could be queried with Astroquery.)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "url = 'https://archive.stsci.edu/pub/hlsp/xdf/hlsp_xdf_hst_acswfc-60mas_hudf_f435w_v1_sci.fits'\n", + "with fits.open(url) as hdulist:\n", + " hdulist.info()\n", + " data = hdulist[0].data\n", + " header = hdulist[0].header" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As explained in the [previous notebook](01_photutils_background_estimation.ipynb) on background estimation, it is important to **mask** these data, as a large portion of the values are equal to zero. We will mask out the non-data portions of the image array, so all of those pixels that have a value of zero don't interfere with our statistics and analyses of the data. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define the mask\n", + "mask = data == 0" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Throughout this notebook, we are going to store our images in Python using a `CCDData` object (see [Astropy documentation](http://docs.astropy.org/en/stable/nddata/index.html#ccddata-class-for-images)), which contains a `numpy` array in addition to metadata such as uncertainty, masks, or units. In this case, our data is in electrons (counts) per second." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "unit = u.electron/ u.s\n", + "xdf_image = CCDData(data, unit=unit, meta=header, mask=mask)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's look at the data:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Set up the figure with subplots\n", + "fig, ax1 = plt.subplots(1, 1, figsize=(8, 8))\n", + "\n", + "# Set up the normalization and colormap\n", + "norm_image = ImageNormalize(vmin=1e-4, vmax=5e-2, stretch=LogStretch(), clip=False)\n", + "cmap = plt.get_cmap('viridis')\n", + "cmap.set_over(cmap.colors[-1])\n", + "cmap.set_under(cmap.colors[0])\n", + "cmap.set_bad('white') # Show masked data as white\n", + "xdf_image_clipped = np.clip(xdf_image, 1e-4, None) # clip to plot with logarithmic stretch\n", + "\n", + "# Plot the data\n", + "fitsplot = ax1.imshow(np.ma.masked_where(xdf_image.mask, xdf_image_clipped), \n", + " norm=norm_image, cmap=cmap)\n", + "\n", + "# Define the colorbar and fix the labels\n", + "cbar = plt.colorbar(fitsplot, fraction=0.046, pad=0.04, ticks=LogLocator(subs=range(10)))\n", + "labels = ['$10^{-4}$'] + [''] * 8 + ['$10^{-3}$'] + [''] * 8 + ['$10^{-2}$']\n", + "cbar.ax.set_yticklabels(labels)\n", + "\n", + "# Define labels\n", + "cbar.set_label(r'Flux Count Rate ({})'.format(xdf_image.unit.to_string('latex')), \n", + " rotation=270, labelpad=30)\n", + "ax1.set_xlabel('X (pixels)')\n", + "ax1.set_ylabel('Y (pixels)')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*Tip: Double-click on any inline plot to zoom in.*" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## Source Detection with `DAOStarFinder`" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "With the `DAOStarFinder` [class](http://photutils.readthedocs.io/en/stable/api/photutils.DAOStarFinder.html), `photutils` provides users with an easy application of the popular [DAOFIND](http://stsdas.stsci.edu/cgi-bin/gethelp.cgi?daofind) algorithm ([Stetson 1987, PASP 99, 191](http://adsabs.harvard.edu/abs/1987PASP...99..191S)), originally developed at the Dominion Astrophysical Observatory. \n", + "\n", + "This algorithm detects sources by:\n", + "* Searching for local maxima\n", + "* Selecting only sources with peak amplitude above a defined threshold\n", + "* Selecting sources with sizes and shapes that match a 2-D Gaussian kernel (circular or elliptical)\n", + "\n", + "It returns:\n", + "* Location of the source centroid\n", + "* Parameters reflecting the source's sharpness and roundness\n", + "\n", + "Generally, the threshold that source detection algorithms use is defined as a multiple of the standard deviation. So first, we need to calculate statistics for the data:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mean, median, std = sigma_clipped_stats(xdf_image.data, sigma=3.0, maxiters=5, mask=xdf_image.mask)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now, let's run the `DAOStarFinder` algorithm on our data and see what it finds. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from photutils import DAOStarFinder" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "daofind = DAOStarFinder(fwhm=5.0, threshold=20.*std)\n", + "sources_dao = daofind(xdf_image * ~xdf_image.mask) \n", + "print(sources_dao)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Set up the figure with subplots\n", + "fig, ax1 = plt.subplots(1, 1, figsize=(8, 8))\n", + "\n", + "# Plot the data\n", + "fitsplot = ax1.imshow(np.ma.masked_where(xdf_image.mask, xdf_image_clipped), \n", + " norm=norm_image, cmap=cmap)\n", + "ax1.scatter(sources_dao['xcentroid'], sources_dao['ycentroid'], s=30, marker='o', \n", + " lw=1, alpha=0.7, facecolor='None', edgecolor='r')\n", + "\n", + "# Define the colorbar and fix the labels\n", + "cbar = plt.colorbar(fitsplot, fraction=0.046, pad=0.04, ticks=LogLocator(subs=range(10)))\n", + "labels = ['$10^{-4}$'] + [''] * 8 + ['$10^{-3}$'] + [''] * 8 + ['$10^{-2}$']\n", + "cbar.ax.set_yticklabels(labels)\n", + "\n", + "# Define labels\n", + "cbar.set_label(r'Flux Count Rate ({})'.format(xdf_image.unit.to_string('latex')), \n", + " rotation=270, labelpad=30)\n", + "ax1.set_xlabel('X (pixels)')\n", + "ax1.set_ylabel('Y (pixels)')\n", + "ax1.set_title('DAOFind Sources')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's randomly pull out some of these sources to get a closer look at them." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, axs = plt.subplots(3,3, figsize=(3, 3))\n", + "plt.subplots_adjust(wspace=0.1, hspace=0.1)\n", + "\n", + "cutout_size = 20\n", + "\n", + "srcs = np.random.permutation(sources_dao)[:axs.size]\n", + "for ax, src in zip(axs.ravel(), srcs):\n", + " slc = (slice(int(src['ycentroid'] - cutout_size), int(src['ycentroid'] + cutout_size)),\n", + " slice(int(src['xcentroid'] - cutout_size), int(src['xcentroid'] + cutout_size)))\n", + " ax.imshow(xdf_image_clipped[slc], norm=norm_image)\n", + " ax.text(2, 2, str(src['id']), color='w', va='top')\n", + " ax.set_xticks([])\n", + " ax.set_yticks([])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Exercises:


\n", + "\n", + "Re-run the `DAOStarFinder` algorithm with a smaller threshold (like 5σ), and plot the sources that it finds. Do the same, but with a larger threshold (like 100σ). How did changing the threshold affect the results?\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Source Detection with `IRAFStarFinder`" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Similarly to `DAOStarFinder`, `IRAFStarFinder` is a class that implements a pre-existing algorithm that is widely used within the astronomical community. This class uses the `starfind` [algorithm](http://stsdas.stsci.edu/cgi-bin/gethelp.cgi?starfind) that was originally part of IRAF.\n", + "\n", + "`IRAFStarFinder` is fundamentally similar to `DAOStarFinder` in that it detects sources by finding local maxima above a certain threshold that match a Gaussian kernel. However, `IRAFStarFinder` differs in the following ways:\n", + "* Does not allow users to specify an elliptical Gaussian kernel\n", + "* Uses image moments to calculate the centroids, roundness, and sharpness of objects\n", + "\n", + "Let's run the `IRAFStarFinder` algorithm on our data, with the same FWHM and threshold, and see how its results differ from `DAOStarFinder`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from photutils import IRAFStarFinder" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "iraffind = IRAFStarFinder(fwhm=5.0, threshold=20.*std)\n", + "sources_iraf = iraffind(xdf_image * ~xdf_image.mask) \n", + "print(sources_iraf)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Set up the figure with subplots\n", + "fig, ax1 = plt.subplots(1, 1, figsize=(8, 8))\n", + "\n", + "# Plot the data\n", + "fitsplot = ax1.imshow(np.ma.masked_where(xdf_image.mask, xdf_image_clipped), \n", + " norm=norm_image, cmap=cmap)\n", + "ax1.scatter(sources_iraf['xcentroid'], sources_iraf['ycentroid'], s=30, marker='o', \n", + " lw=1, alpha=0.7, facecolor='None', edgecolor='r')\n", + "\n", + "# Define the colorbar and fix the labels\n", + "cbar = plt.colorbar(fitsplot, fraction=0.046, pad=0.04, ticks=LogLocator(subs=range(10)))\n", + "labels = ['$10^{-4}$'] + [''] * 8 + ['$10^{-3}$'] + [''] * 8 + ['$10^{-2}$']\n", + "cbar.ax.set_yticklabels(labels)\n", + "\n", + "# Define labels\n", + "cbar.set_label(r'Flux Count Rate ({})'.format(xdf_image.unit.to_string('latex')), \n", + " rotation=270, labelpad=30)\n", + "ax1.set_xlabel('X (pixels)')\n", + "ax1.set_ylabel('Y (pixels)')\n", + "ax1.set_title('IRAFFind Sources')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Again, let's randomly select some sources for a closer look:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, axs = plt.subplots(3,3, figsize=(3, 3))\n", + "plt.subplots_adjust(wspace=0.1, hspace=0.1)\n", + "\n", + "cutout_size = 20\n", + "\n", + "srcs = np.random.permutation(sources_iraf)[:axs.size]\n", + "for ax, src in zip(axs.ravel(), srcs):\n", + " slc = (slice(int(src['ycentroid'] - cutout_size), int(src['ycentroid'] + cutout_size)),\n", + " slice(int(src['xcentroid'] - cutout_size), int(src['xcentroid'] + cutout_size)))\n", + " ax.imshow(xdf_image_clipped[slc], norm=norm_image)\n", + " ax.text(2, 2, str(src['id']), color='w', va='top')\n", + " ax.set_xticks([])\n", + " ax.set_yticks([])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Exercises:


\n", + "\n", + "Re-run the `IRAFStarFinder` algorithm with a smaller full-width-half-max (FWHM) – try 3 pixels – and plot the sources that it finds. Do the same, but with a larger FWHM (like 10 pixels). How did changing the FWHM affect the results? What astronomical objects might be better captures by smaller FWHM? Larger?\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Note: Comparing `DAOStarFinder` and `IRAFStarFinder`" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You might have noticed that the `IRAFStarFinder` algorithm only found 211 sources in our data - 14% of what `DAOStarFinder` found. Why is this?\n", + "\n", + "The answer comes down to the default settings for the two algorithms: (1) there are differences in the upper and lower bounds on the requirements for source roundness and sharpness, and (2) `IRAFStarFinder` includes a minimum separation between sources that `DAOStarFinder` does not have:\n", + "\n", + "| | `IRAFStarFinder` | `DAOStarFinder` |\n", + "|------|------|------|\n", + "| sharplo | 0.5 | 0.2 |\n", + "| sharphi | 2.0 | 1.0 |\n", + "| roundlo | 0.0 | -1.0 |\n", + "| roundhi | 0.2 | 1.0 |\n", + "| minsep_fwhm | 1.5 * FWHM | N/A |\n", + "\n", + "Thinking about this, *it then makes sense* that `IRAFStarFinder` would find fewer sources. It has stricter restrictions on source roundness, meaning that it eliminates more elliptical galactic sources (this is the eXtreme Deep Field, after all!), and the minimum separation requirement further rules out sources that are too close to one another.\n", + "\n", + "If we set all these parameters to be equivalent, though, we should find much better agreement between the two methods:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "iraffind_match = IRAFStarFinder(fwhm=5.0, threshold=20.*std, \n", + " sharplo=0.2, sharphi=1.0, \n", + " roundlo=-1.0, roundhi=1.0,\n", + " minsep_fwhm=0.0)\n", + "sources_iraf_match = iraffind_match(xdf_image * ~xdf_image.mask) \n", + "print(sources_iraf_match)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The number of detected sources are in much better agreement now - 1415 versus 1470 - but the improved agreement can also be seen by plotting the location of these sources:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Set up the figure with subplots\n", + "fig, [ax1, ax2] = plt.subplots(1, 2, figsize=(12, 6))\n", + "plt.tight_layout()\n", + "\n", + "# Plot the DAOStarFinder data\n", + "fitsplot = ax1.imshow(np.ma.masked_where(xdf_image.mask, xdf_image_clipped), \n", + " norm=norm_image, cmap=cmap)\n", + "ax1.scatter(sources_dao['xcentroid'], sources_dao['ycentroid'], s=30, marker='o', \n", + " lw=1, alpha=0.7, facecolor='None', edgecolor='r')\n", + "ax1.set_xlabel('X (pixels)')\n", + "ax1.set_ylabel('Y (pixels)')\n", + "ax1.set_title('DAOStarFinder Sources')\n", + "\n", + "# Plot the IRAFStarFinder data\n", + "fitsplot = ax2.imshow(np.ma.masked_where(xdf_image.mask, xdf_image_clipped), \n", + " norm=norm_image, cmap=cmap)\n", + "ax2.scatter(sources_iraf_match['xcentroid'], sources_iraf_match['ycentroid'], \n", + " s=30, marker='o', lw=1, alpha=0.7, facecolor='None', edgecolor='r')\n", + "ax2.set_xlabel('X (pixels)')\n", + "ax2.set_title('IRAFStarFinder Sources')\n", + "\n", + "# Define the colorbar\n", + "cbar_ax = fig.add_axes([1, 0.09, 0.03, 0.87])\n", + "cbar = plt.colorbar(fitsplot, cbar_ax, ticks=LogLocator(subs=range(10)))\n", + "labels = ['$10^{-4}$'] + [''] * 8 + ['$10^{-3}$'] + [''] * 8 + ['$10^{-2}$']\n", + "cbar.ax.set_yticklabels(labels)\n", + "cbar.set_label(r'Flux Count Rate ({})'.format(xdf_image.unit.to_string('latex')), \n", + " rotation=270, labelpad=30)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Take this example as reminder to be mindful when selecting a source detection algorithm, and when defining algorithm parameters! Don't be afraid to play around with the parameters and investigate how that affects your results." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Source Detection with `find_peaks`" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For more general source detection cases that do not require comparison with models, `photutils` offers the `find_peaks` function. \n", + "\n", + "This function simply finds sources by identifying local maxima above a given threshold and separated by a given distance, rather than trying to fit data to a given model. Unlike the previous detection algorithms, `find_peaks` does not necessarily calculate objects' centroids. Unless the `subpixel` argument is set to `True`, `find_peaks` will return just the integer value of the peak pixel for each source.\n", + "\n", + "This algorithm is particularly useful for identifying non-stellar sources or heavily distorted sources in image data.\n", + "\n", + "Let's see how it does:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from photutils import find_peaks" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sources_findpeaks = find_peaks(xdf_image.data, mask=xdf_image.mask, \n", + " threshold=20.*std, box_size=30, subpixel=True) \n", + "print(sources_findpeaks)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Set up the figure with subplots\n", + "fig, ax1 = plt.subplots(1, 1, figsize=(8, 8))\n", + "\n", + "# Plot the data\n", + "fitsplot = ax1.imshow(np.ma.masked_where(xdf_image.mask, xdf_image_clipped), \n", + " norm=norm_image, cmap=cmap)\n", + "ax1.scatter(sources_findpeaks['x_peak'], sources_findpeaks['y_peak'], s=30, marker='o', \n", + " lw=1, alpha=0.7, facecolor='None', edgecolor='r')\n", + "\n", + "# Define the colorbar\n", + "cbar = plt.colorbar(fitsplot, fraction=0.046, pad=0.04, ticks=LogLocator(subs=range(10)))\n", + "labels = ['$10^{-4}$'] + [''] * 8 + ['$10^{-3}$'] + [''] * 8 + ['$10^{-2}$']\n", + "cbar.ax.set_yticklabels(labels)\n", + "\n", + "# Define labels\n", + "cbar.set_label(r'Flux Count Rate ({})'.format(xdf_image.unit.to_string('latex')), \n", + " rotation=270, labelpad=30)\n", + "ax1.set_xlabel('X (pixels)')\n", + "ax1.set_ylabel('Y (pixels)')\n", + "ax1.set_title('find\\_peaks Sources')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And a closer look:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, axs = plt.subplots(3,3, figsize=(3, 3))\n", + "plt.subplots_adjust(wspace=0.1, hspace=0.1)\n", + "\n", + "cutout_size = 20\n", + "\n", + "srcs = np.random.permutation(sources_findpeaks)[:axs.size]\n", + "for ax, src in zip(axs.ravel(), srcs):\n", + " slc = (slice(int(src['y_peak'] - cutout_size), int(src['y_peak'] + cutout_size)),\n", + " slice(int(src['x_peak'] - cutout_size), int(src['x_peak'] + cutout_size)))\n", + " ax.imshow(xdf_image_clipped[slc], norm=norm_image)\n", + " src_id = np.where((sources_findpeaks['x_peak'] == src['x_peak']) & \n", + " (sources_findpeaks['y_peak'] == src['y_peak']))[0][0]\n", + " ax.text(2, 2, str(src_id), color='w', va='top')\n", + " ax.set_xticks([])\n", + " ax.set_yticks([])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Comparing Detection Methods\n", + "\n", + "Let's compare how each of these different strategies did.\n", + "\n", + "First, how many sources did each method find?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print('''DAOStarFinder: {} sources\n", + "IRAFStarFinder: {} sources\n", + "find_peaks: {} sources'''.format(len(sources_dao), len(sources_iraf), len(sources_findpeaks)))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, how many of these sources match? We can answer this question by using [sets](https://docs.python.org/3/library/stdtypes.html#set-types-set-frozenset) to compare the centroids of the different sources (rounding to the first decimal place)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Make lists of centroid coordinates\n", + "centroids_dao = [(x, y) for x, y in sources_dao['xcentroid', 'ycentroid']]\n", + "centroids_iraf = [(x, y) for x, y in sources_iraf['xcentroid', 'ycentroid']]\n", + "centroids_findpeaks = [(x, y) for x, y in sources_findpeaks['x_centroid', 'y_centroid']]\n", + "\n", + "# Round those coordinates to the first decimal place and convert them to be sets\n", + "rounded_centroids_dao = set([(round(x, 1), round(y, 1)) for x, y in centroids_dao])\n", + "rounded_centroids_iraf = set([(round(x, 1), round(y, 1)) for x, y in centroids_iraf])\n", + "rounded_centroids_findpeaks = set([(round(x, 1), round(y, 1)) for x, y in centroids_findpeaks])\n", + "\n", + "# Examine the intersections of different sets to determine which sources are shared\n", + "all_match = rounded_centroids_dao.intersection(rounded_centroids_iraf).intersection(rounded_centroids_findpeaks)\n", + "dao_iraf_match = rounded_centroids_dao.intersection(rounded_centroids_iraf)\n", + "dao_findpeaks_match = rounded_centroids_dao.intersection(rounded_centroids_findpeaks)\n", + "iraf_findpeaks_match = rounded_centroids_iraf.intersection(rounded_centroids_findpeaks)\n", + "\n", + "print('''Matching sources found by:\n", + " All methods: {}\n", + " DAOStarFinder & IRAFStarFinder: {}\n", + " DAOStarFinder & find_peaks: {}\n", + " IRAFStarFinder & find_peaks: {}'''\n", + " .format(len(all_match), len(dao_iraf_match), \n", + " len(dao_findpeaks_match), len(iraf_findpeaks_match)))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And just for fun, let's plot these matching sources. (The colors chosen to represent different sets are from [Paul Tol's guide for accessible color schemes](https://personal.sron.nl/~pault/).)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Set up the figure with subplots\n", + "fig, ax1 = plt.subplots(1, 1, figsize=(8, 8))\n", + "\n", + "# Plot the data\n", + "fitsplot = ax1.imshow(np.ma.masked_where(xdf_image.mask, xdf_image_clipped), norm=norm_image)\n", + "ax1.scatter([x for x, y in list(dao_findpeaks_match)], [y for x, y in list(dao_findpeaks_match)],\n", + " s=30, marker='s', lw=1.5, facecolor='None', edgecolor='#EE7733',\n", + " label='Found by DAO \\& find\\_peaks')\n", + "ax1.scatter([x for x, y in list(dao_iraf_match)], [y for x, y in list(dao_iraf_match)],\n", + " s=30, marker='D', lw=1.5, facecolor='None', edgecolor='#EE3377',\n", + " label='Found by DAO \\& IRAF')\n", + "ax1.scatter([x for x, y in list(iraf_findpeaks_match)], [y for x, y in list(iraf_findpeaks_match)],\n", + " s=30, marker='o', lw=1.5, facecolor='None', edgecolor='#0077BB',\n", + " label='Found by IRAF \\& find\\_peaks')\n", + "ax1.scatter([x for x, y in list(all_match)], [y for x, y in list(all_match)],\n", + " s=30, marker='o', lw=1.2, linestyle=':',facecolor='None', edgecolor='#BBBBBB',\n", + " label='Found by all methods')\n", + "ax1.legend()\n", + "\n", + "# Define the colorbar\n", + "cbar = plt.colorbar(fitsplot, fraction=0.046, pad=0.04, ticks=LogLocator(subs=range(10)))\n", + "labels = ['$10^{-4}$'] + [''] * 8 + ['$10^{-3}$'] + [''] * 8 + ['$10^{-2}$']\n", + "cbar.ax.set_yticklabels(labels)\n", + "\n", + "# Define labels\n", + "cbar.set_label(r'Flux Count Rate ({})'.format(xdf_image.unit.to_string('latex')), \n", + " rotation=270, labelpad=30)\n", + "ax1.set_xlabel('X (pixels)')\n", + "ax1.set_ylabel('Y (pixels)')\n", + "ax1.set_title('Sources Found by Different Methods')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*Remember that you can double-click on the plot to zoom in and look around!*" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Custom Detection Algorithms\n", + "\n", + "If none of the algorithms we've reviewed above do exactly what you need, `photutils` also provides infrastructure for you to generate and use your own source detection algorithm: the `StarFinderBase` object can be inherited and used to develop new star-finding classes. Take a look at the [documentation](https://photutils.readthedocs.io/en/latest/api/photutils.detection.StarFinderBase.html#photutils.detection.StarFinderBase) for more information.\n", + "\n", + "If you do go that route, remember that `photutils` is open-developed; you would be very welcome to [open a pull request](https://github.com/astropy/photutils/blob/master/CONTRIBUTING.rst) and incorporate your new star finder into the `photutils` source code - for everyone to use!" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "\n", + "## Image Segmentation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Beyond traditional source detection methods, an additional option for identifying sources in image data is a process called **image segmentation**. This method identifies and labels contiguous (connected) objects within an image. \n", + "\n", + "You might have noticed that, in the previous source detection algorithms, large and extended sources are often incorrectly identified as more than one source. Segmentation would label all the pixels within a large galaxy as belonging to the same object, and would allow the user to then measure the photometry, centroid, and morphology of the entire object at once.\n", + "\n", + "#### Creating a `SegmentationImage`\n", + "\n", + "In `photutils`, image segmentation maps are created using the threshold method in the `detect_sources()` function. This method identifies all of the objects in the data that have signals above a determined **`threshold`** (usually defined as a multiple of the standard deviation) and that have more than a defined number of adjoining pixels, **`npixels`**. The data can also optionally be smoothed using a kernel, **`filter_kernel`**, before applying the threshold cut.\n", + "\n", + "The `detect_sources()` function returns a `SegmentationImage` object: an array in which each object is labeled with an integer. As a simple example, a segmentation map containing two distinct sources might look like this:\n", + "\n", + "```\n", + "0 0 0 0 0 0 0 0 0 0\n", + "0 1 1 0 0 0 0 0 0 0\n", + "1 1 1 1 1 0 0 0 2 0\n", + "1 1 1 1 0 0 0 2 2 2\n", + "1 1 1 0 0 0 2 2 2 2\n", + "1 1 1 1 0 0 0 2 2 0\n", + "1 1 0 0 0 0 2 2 0 0\n", + "0 1 0 0 0 0 2 0 0 0\n", + "0 0 0 0 0 0 0 0 0 0\n", + "```\n", + "where all of the pixels labeled `1` belong to the first source, all those labeled `2` belong to the second, and all null pixels are designated to be background.\n", + "\n", + "Let's see what the segmentation map for our XDF data will look like." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from photutils import detect_sources\n", + "from photutils.utils import random_cmap" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define threshold and minimum object size\n", + "threshold = 5. * std\n", + "npixels = 15\n", + "\n", + "# Create a segmentation image\n", + "segm = detect_sources(xdf_image.data, threshold, npixels)\n", + "\n", + "print('Found {} sources'.format(segm.max_label))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Set up the figure with subplots\n", + "fig, [ax1, ax2] = plt.subplots(1, 2, figsize=(12, 6))\n", + "plt.tight_layout()\n", + "\n", + "# Plot the original data\n", + "fitsplot = ax1.imshow(np.ma.masked_where(xdf_image.mask, xdf_image_clipped), \n", + " norm=norm_image, cmap=cmap)\n", + "ax1.set_ylabel('Y (pixels)')\n", + "ax1.set_title('Original Data')\n", + "\n", + "# Plot the segmentation image\n", + "rand_cmap = random_cmap(random_state=12345)\n", + "rand_cmap.set_under('black')\n", + "segplot = ax2.imshow(np.ma.masked_where(xdf_image.mask, segm), vmin=1, cmap=rand_cmap)\n", + "ax2.set_xlabel('X (pixels)')\n", + "ax2.set_title('Segmentation Map')\n", + "\n", + "# Define the colorbar\n", + "cbar_ax = fig.add_axes([1, 0.09, 0.03, 0.87])\n", + "cbar = plt.colorbar(segplot, cbar_ax)\n", + "cbar.set_label('Object Label', rotation=270, labelpad=30)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Compare the sources in original data to those in the segmentation image. Each color in the segmentation map denotes a separate source.\n", + "\n", + "You can easily see that larger galaxies are shown in the segmentation map as contiguous objects of the same color - for example, the two yellow and pink galaxies near (1200, 2500). Each pixel containing light from the same galaxy has been labeled as belonging to the same object." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For a closer look, let's see what the sources we found with `find_peaks` look like in this segmentation map:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "fig, axs = plt.subplots(3,3, figsize=(3, 3))\n", + "plt.subplots_adjust(wspace=0.1, hspace=0.1)\n", + "\n", + "cutout_size = 20\n", + "\n", + "srcs = np.random.permutation(sources_findpeaks)[:axs.size]\n", + "for ax, src in zip(axs.ravel(), srcs):\n", + " slc = (slice(int(src['y_peak'] - cutout_size), int(src['y_peak'] + cutout_size)),\n", + " slice(int(src['x_peak'] - cutout_size), int(src['x_peak'] + cutout_size)))\n", + " ax.imshow(segm.data[slc], cmap=rand_cmap, vmin=1, vmax=len(sources_findpeaks))\n", + " src_id = np.where((sources_findpeaks['x_peak'] == src['x_peak']) & \n", + " (sources_findpeaks['y_peak'] == src['y_peak']))[0][0]\n", + " ax.text(2, 2, str(src_id), color='w', va='top')\n", + " ax.set_xticks([])\n", + " ax.set_yticks([])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Exercises:


\n", + "\n", + "Recompute the `SegmentationImage`, but alter the threshold and the minimum number of pixels in a source. How does changing the threshold affect the results? What about changing the number of pixels?\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Analyzing `source_properties`\n", + "\n", + "Once we have a `SegmentationImage` object, `photutils` provides many powerful tools to manipulate and analyze the identified objects. \n", + "\n", + "Individual objects within the segmentation map can be altered using methods such as `relabel` to change the labels of objects, `remove_labels` to remove objects, or `deblend_sources` to separating overlapping sources that were incorrectly labeled as one source.\n", + "\n", + "However, perhaps the most powerful aspect of the `SegmentationImage` is the ability to create a catalog using `source_properties` to measure the centroids, photometry, and morphology of the detected objects:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from photutils import source_properties" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "catalog = source_properties(xdf_image.data, segm)\n", + "table = catalog.to_table()\n", + "print(table)\n", + "print(table.colnames)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Creating apertures from segmentation data\n", + "\n", + "We can use this information to create isophotal ellipses for each identified source. These ellipses can also later be used as photometric apertures!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from photutils import EllipticalAperture" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Define the approximate isophotal extent\n", + "r = 4.\n", + "\n", + "# Create the apertures\n", + "apertures = []\n", + "for obj in catalog:\n", + " position = (obj.xcentroid.value, obj.ycentroid.value)\n", + " a = obj.semimajor_axis_sigma.value * r\n", + " b = obj.semiminor_axis_sigma.value * r\n", + " theta = obj.orientation.value\n", + " apertures.append(EllipticalAperture(position, a, b, theta=theta))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Set up the figure with subplots\n", + "fig, ax1 = plt.subplots(1, 1, figsize=(8, 8))\n", + "\n", + "# Plot the data\n", + "fitsplot = ax1.imshow(np.ma.masked_where(xdf_image.mask, xdf_image_clipped), \n", + " norm=norm_image, cmap=cmap)\n", + "\n", + "# Plot the apertures\n", + "for aperture in apertures:\n", + " aperture.plot(color='red', lw=1, alpha=0.7, ax=ax1)\n", + "\n", + "# Define the colorbar\n", + "cbar = plt.colorbar(fitsplot, fraction=0.046, pad=0.04, ticks=LogLocator(subs=range(10)))\n", + "labels = ['$10^{-4}$'] + [''] * 8 + ['$10^{-3}$'] + [''] * 8 + ['$10^{-2}$']\n", + "cbar.ax.set_yticklabels(labels)\n", + "\n", + "# Define labels\n", + "cbar.set_label(r'Flux Count Rate ({})'.format(xdf_image.unit.to_string('latex')), \n", + " rotation=270, labelpad=30)\n", + "ax1.set_xlabel('X (pixels)')\n", + "ax1.set_ylabel('Y (pixels)')\n", + "ax1.set_title('Segmentation Image Apertures')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Exercises:


\n", + "\n", + "Play with the isophotal extent of the elliptical apertures (defined above as `r`). Observe how changing this value affects the apertures that are created.\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "It is clear that using `photutils` for image segmentation can allow users to generate highly customized apertures - great for complex data that contain many different kinds of celestial sources." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## Conclusions\n", + "\n", + "The `photutils` package provides users with a variety of methods for detecting sources in their data, from familar algorithms such as `DAOFind` and `starfind`, to more complex and customizable image segmentation algorithms. These methods allow for easy creation of a diverse array of apertures that can be used for photometric analysis.\n", + "\n", + "**To continue with this `photutils` tutorial, go on to the [aperture photometry](03_photutils_aperture_photometry.ipynb) or [PSF photometry notebook](04_photutils_psf_photometry.ipynb).**" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## Additional Resources\n", + "For more examples and details, please visit the [photutils](http://photutils.readthedocs.io/en/stable/index.html) documentation." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "## About this Notebook\n", + "**Authors:** Lauren Chambers (lchambers@stsci.edu), Erik Tollerud (etollerud@stsci.edu), Tom Wilson (towilson@stsci.edu)\n", + "
**Updated:** May 2019" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "[Top of Page](#title_ID)\n", + "\"STScI" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.6.8" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/requirements.txt b/requirements.txt index 6b03f2e..367e4e4 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,7 @@ astropy +jupyter +photutils +matplotlib nbconvert nbformat -jupyter +numpy