{ "cells": [ { "cell_type": "markdown", "id": "9699e83a", "metadata": {}, "source": [ "# Quick Start\n", "\n", "This section provides an introduction to using `PetroFit`. \n", "It also highlights essential functions in conjunction with `Photutils` and `Astropy`. \n", "For in-depth explanations, refer to the respective sections of interest within this documentation.\n", "\n", "To start with `PetroFit`, simply import it as follows:" ] }, { "cell_type": "code", "execution_count": null, "id": "bc6d9298", "metadata": {}, "outputs": [], "source": [ "import petrofit as pf" ] }, { "cell_type": "markdown", "id": "a68cbd8e", "metadata": {}, "source": [ "## Loading Example Data\n", "\n", "The following data is a cutout of a group of bright galaxies in Abell 2744 (located at `(3.596248, -30.388517)`). The original data was acquired by the [Hubble Frontier Fields](https://frontierfields.org) team via the WFC3 instrument in the `F105W` filter and can be directly downloaded from the [Mikulski Archive for Space Telescopes](https://archive.stsci.edu/pub/hlsp/frontier/abell2744/images/hst/v1.0/hlsp_frontier_hst_wfc3-60mas_abell2744_f105w_v1.0_drz.fits). The cutout image used in this documentation can be found in the git repository at the following path `petrofit/docs/data/abell_2744_dwarf_galaxy_f105w.fits.gz`.\n", "\n", "### Loading Image \n", "You can use `astropy`'s ``CCDData`` to load the example data and visualize it through `matplotlib`. Note that ``CCDData.read`` does not return a 2D array, but rather a ``CCDData`` instance which contains the image array, header, and WCS. To access the image array stored in the ``CCDData`` use the ``data`` attribute (i.e. ``CCDData.data`` as shown in the ``plt.imshow`` command below)." ] }, { "cell_type": "code", "execution_count": null, "id": "33eab6a6", "metadata": {}, "outputs": [], "source": [ "from astropy.nddata import CCDData\n", "\n", "image = CCDData.read('data/abell_2744_dwarf_galaxy_f105w.fits.gz')" ] }, { "cell_type": "code", "execution_count": null, "id": "7c75fd46", "metadata": { "nbsphinx": "hidden" }, "outputs": [], "source": [ "# Hidden cell\n", "\n", "%matplotlib inline\n", "\n", "# Stop Fit Model to Data section warnings\n", "import warnings\n", "warnings.filterwarnings('ignore', append=True)" ] }, { "cell_type": "code", "execution_count": null, "id": "f2f045f6", "metadata": { "scrolled": false }, "outputs": [], "source": [ "import numpy as np\n", "from matplotlib import pyplot as plt\n", "\n", "plt.rcParams['figure.figsize'] = [6, 6]\n", "plt.rcParams['image.origin'] = 'lower'\n", "plt.rcParams['font.size'] = 12\n", "\n", "vmax = 0.05 # vmax for matplotlib imshow\n", "vmin = - vmax \n", "\n", "plt.imshow(image.data, vmin=vmin, vmax=vmax)\n", "plt.title(\"Abell 2744 Galaxies\")\n", "plt.xlabel(\"Pixels\")\n", "plt.ylabel(\"Pixels\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "1b5204cf", "metadata": {}, "source": [ "### Loading RMS Image\n", "Since we only want the rms image array, we use `astropy`'s ``io.fits.getdata`` function as follows:" ] }, { "cell_type": "code", "execution_count": null, "id": "61e43267", "metadata": {}, "outputs": [], "source": [ "from astropy.io import fits \n", "rms = fits.getdata('data/abell_2744_dwarf_galaxy_f105w_rms.fits.gz')" ] }, { "cell_type": "code", "execution_count": null, "id": "6a9f164d", "metadata": {}, "outputs": [], "source": [ "plt.imshow(rms)\n", "plt.title(\"RMS Image\")\n", "plt.xlabel(\"Pixels\")\n", "plt.ylabel(\"Pixels\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "004686be", "metadata": {}, "source": [ "### Making Cutouts\n", "\n", "You can use `astropy`'s `Cutout2D` function to make cutouts of sources. To access the data (image array), use the ``data`` attribute (i.e. ``Cutout2D.data`` as shown in the ``plt.imshow`` command below). Note that ``position`` can be a `SkyCoord` if you provide the `Cutout2D` function a WCS." ] }, { "cell_type": "code", "execution_count": null, "id": "fd575cb6", "metadata": {}, "outputs": [], "source": [ "from astropy.nddata import Cutout2D\n", "\n", "# Make cutout image, centerd at (100, 100) pixels, 40 pixels in size\n", "cutout_image = Cutout2D(image, position=(100,100), size=40)\n", "\n", "# Make cutout rms, centerd at (100, 100) pixels, 40 pixels in size\n", "cutout_rms = Cutout2D(rms, position=(100,100), size=40)\n", "\n", "# Plot cutouts\n", "# ------------\n", "fig, axs = plt.subplots(1,2, figsize=(12, 6))\n", "\n", "plt.sca(axs[0])\n", "plt.imshow(cutout_image.data, vmin=vmin, vmax=vmax)\n", "plt.title(\"Cutout Galaxy\")\n", "plt.xlabel(\"Pixels\")\n", "plt.ylabel(\"Pixels\")\n", "\n", "plt.sca(axs[1])\n", "plt.imshow(cutout_rms.data)\n", "plt.title(\"Cutout RMS\")\n", "plt.xlabel(\"Pixels\")\n", "plt.ylabel(\"Pixels\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "0cf2d635", "metadata": {}, "source": [ "## Defining AstroPy Sérsic Models" ] }, { "cell_type": "markdown", "id": "7af930fe", "metadata": {}, "source": [ "This example shows how to define a 2D Sérsic model using `astropy`. We fill in our initial guess for the parameters (or correct parameters if we know them) when initializing the ``Sersic2D`` object. \n", "\n", "To assist users, PetroFit offers a function named `get_default_sersic_bounds` that provides a default set of bounds for the Sérsic model parameters. This function is particularly useful when you're fitting data as these bounds help in constraining the fitting parameter space. \n", "\n", "The dictionary returned by `get_default_sersic_bounds` has parameter names as its keys. The values are tuples where the first element indicates the minimum allowed value and the second element indicates the maximum. In situations where there's no constraint on a parameter's value, `None` is used.\n", "\n", "For instance, the dictionary looks like this:\n", "\n", "\n", "```python\n", "# get_default_sersic_bounds returns:\n", "bounds = {\n", " 'amplitude': (0., None),\n", " 'r_eff': (0, None),\n", " 'n': (0, 10),\n", " 'ellip': (0, 1),\n", " 'theta': (-2 * np.pi, 2 * np.pi),\n", "}\n", "```" ] }, { "cell_type": "code", "execution_count": null, "id": "0b50d7c6", "metadata": {}, "outputs": [], "source": [ "from astropy.modeling import models\n", "\n", "sersic_model = models.Sersic2D(\n", "\n", " amplitude=10, # Intensity at r_eff\n", " r_eff=1, # Effective or half-lilght radius\n", " n=4, # Sersic index\n", " x_0=20, # center of model in the x direction\n", " y_0=20, # center of model in the y direction\n", " ellip=0.1, # Ellipticity\n", " theta=0.0, # Rotation angle in radians, counterclockwise from the positive x-axis.\n", " \n", " bounds=pf.get_default_sersic_bounds(), # PetroFit parameter bounds\n", ")" ] }, { "cell_type": "markdown", "id": "8ec21f0f", "metadata": {}, "source": [ "To add `x_0` and `y_0` bounds to the default bounds, you can update the dictionary as you would a regular Python dictionary. You can add or update bounds by passing a Python dictionary to the `get_default_sersic_bounds` function as follows:" ] }, { "cell_type": "code", "execution_count": null, "id": "fa21fd37", "metadata": {}, "outputs": [], "source": [ "bound_dict = pf.get_default_sersic_bounds({\n", " 'x_0': (10, 30), \n", " 'y_0': (10, 30)})\n", "bound_dict" ] }, { "cell_type": "markdown", "id": "80ec53bd", "metadata": {}, "source": [ "You can also directly update the model bounds as follows:" ] }, { "cell_type": "code", "execution_count": null, "id": "98ab2630", "metadata": {}, "outputs": [], "source": [ "sersic_model.bounds.update({'x_0': (10, 30), 'y_0': (10, 30)})" ] }, { "cell_type": "markdown", "id": "0d9fee69", "metadata": {}, "source": [ "## Making Compound Models (Combining Models)\n", "\n", "You can combine multiple models to form a compound model by adding, subtracting, multiplying, and dividing individual models. For example, we add the Sérsic model from the last section to itself to form a two-component Sérsic model (notice that the number of parameters double):" ] }, { "cell_type": "code", "execution_count": null, "id": "e710d764", "metadata": {}, "outputs": [], "source": [ "compound_sersic_model = sersic_model + sersic_model" ] }, { "cell_type": "code", "execution_count": null, "id": "1836103d", "metadata": {}, "outputs": [], "source": [ "pf.print_model_params(compound_sersic_model)" ] }, { "cell_type": "markdown", "id": "f1fbe442", "metadata": {}, "source": [ "## Making a PSF Convolved Model\n", "\n", "`PSFConvolvedModel2D` is an extension of `Fittable2DModel` designed to incorporate PSF convolution and facilitate image sampling for models native to `astropy`. This model functions by:\n", "\n", "1. Generating an image representation of the base model.\n", "2. Sampling this image onto a defined grid.\n", "3. If supplied, convolving the model image with a PSF.\n", "\n", "Because `PSFConvolvedModel2D` retains the properties of a `Fittable2DModel`, it can be used for fitting model-generated images against actual data. \n", "As an illustration, this documentation will demonstrate the usage of `PSFConvolvedModel2D` with `astropy`'s `Sersic2D` model we defined above.\n" ] }, { "cell_type": "code", "execution_count": null, "id": "a1d501d2", "metadata": {}, "outputs": [], "source": [ "# Load PSF image (2D array)\n", "PSF = fits.getdata('data/f105w_psf.fits.gz')\n", "\n", "# Normalize PSF \n", "PSF = PSF / PSF.sum()\n", "\n", "# Note that the PSF shape is odd on all sides\n", "print(\"PSF Shape = {}\".format(PSF.shape))\n", "\n", "# Plot PSF and use vmax and vmin to show difraction spikes\n", "plt.imshow(PSF, vmin=0, vmax=PSF.std()/10)\n", "plt.title(\"PSF Image\")\n", "plt.xlabel(\"Pixels\")\n", "plt.ylabel(\"Pixels\")\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "0ff9fa9f", "metadata": {}, "outputs": [], "source": [ "# PSFConvolvedModel2D \n", "psf_sersic_model = pf.PSFConvolvedModel2D(sersic_model, psf=PSF, oversample=4)" ] }, { "cell_type": "markdown", "id": "a308c0df", "metadata": {}, "source": [ "## Converting Models to Images\n", "\n", "To convert any 2D model (Astropy or PetroFit) to an image use the `model_to_image` function." ] }, { "cell_type": "code", "execution_count": null, "id": "fbf6b956", "metadata": {}, "outputs": [], "source": [ "# Size of model image\n", "size = 40\n", "\n", "# sersic model image\n", "model_image = pf.model_to_image(model=sersic_model, size=size)\n", "\n", "# PSF convolved model image \n", "psf_model_image = pf.model_to_image(model=psf_sersic_model, size=size)" ] }, { "cell_type": "markdown", "id": "39f102fe", "metadata": {}, "source": [ "Plot model image" ] }, { "cell_type": "code", "execution_count": null, "id": "6ca12c73", "metadata": { "scrolled": false }, "outputs": [], "source": [ "fig, axs = plt.subplots(1,2, figsize=(12, 6))\n", "\n", "plt.sca(axs[0])\n", "plt.imshow(model_image, vmin=vmin, vmax=vmax)\n", "plt.title('Sérsic Model Image (n=4)')\n", "plt.xlabel(\"Pixels\")\n", "plt.ylabel(\"Pixels\")\n", "\n", "plt.sca(axs[1])\n", "plt.imshow(psf_model_image, vmin=vmin, vmax=vmax)\n", "plt.title('PSF Convolved Sérsic Model Image (n=4)')\n", "plt.xlabel(\"Pixels\")\n", "plt.ylabel(\"Pixels\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "5777eea6", "metadata": {}, "source": [ "## Fitting Model to Image\n", "\n", "We first define a `PSFConvolvedModel2D` model with initial guesses as follows:" ] }, { "cell_type": "code", "execution_count": null, "id": "462d0fd6", "metadata": {}, "outputs": [], "source": [ "sersic_model = models.Sersic2D(\n", " amplitude=0.1, # Intensity at r_eff\n", " r_eff=10, # Effective or half-lilght radius\n", " n=1, # Sersic index\n", " x_0=20, # center of model in the x direction\n", " y_0=20, # center of model in the y direction\n", " ellip=0.1, # Ellipticity\n", " theta=0.0, # Rotation angle in radians, counterclockwise from the positive x-axis.\n", " \n", " bounds=pf.get_default_sersic_bounds({'x_0': (10, 30), 'y_0': (10, 30)}), # Parameter bounds\n", ")\n", "\n", "psf_sersic_model = pf.PSFConvolvedModel2D(sersic_model, psf=PSF, oversample=4)" ] }, { "cell_type": "markdown", "id": "819ecbba", "metadata": {}, "source": [ "Before we fit the image, we compute the weights of each pixel using rms data as follows (please note that weights are optional and set to `None` by default):" ] }, { "cell_type": "code", "execution_count": null, "id": "3a3a82c6", "metadata": {}, "outputs": [], "source": [ "fitting_weights = 1 / cutout_rms.data " ] }, { "cell_type": "markdown", "id": "6ed84fb9", "metadata": {}, "source": [ "Use the `fit_model` function to fit 2D models to images as follows:" ] }, { "cell_type": "code", "execution_count": null, "id": "1f7261ad", "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "fitted_model, fitter = pf.fit_model(\n", " image=cutout_image.data, \n", " model=psf_sersic_model,\n", " weights=fitting_weights, \n", " maxiter=10000,\n", " epsilon=1.4901161193847656e-08,\n", " acc=1e-09,\n", ")" ] }, { "cell_type": "markdown", "id": "9ad84ae6", "metadata": {}, "source": [ "Convert the fitted model into an image" ] }, { "cell_type": "code", "execution_count": null, "id": "f64c5cf0", "metadata": {}, "outputs": [], "source": [ "axs, cbar, model_image, residual_image = pf.plot_fit(fitted_model, cutout_image, \n", " vmax=vmax, vmin=vmin, figsize=[6*3, 6])\n", " \n", "plt.show()" ] }, { "cell_type": "markdown", "id": "92619d47", "metadata": {}, "source": [ "## Fitting Multiple Sources\n", "\n", "If the locations of the sources are known, we can fit all sources at the same time by creating a compound model. Note that `x_0` and `y_0` are known beforehand using photometric centroids. Below, a compound model of 3 Sérsic components is defined and the original image is fit (i.e. not the cutout we have been working with).\n" ] }, { "cell_type": "code", "execution_count": null, "id": "fdd363e4", "metadata": {}, "outputs": [], "source": [ "# Center elliptical galaxy we have been fitting:\n", "galaxy_model_1 = models.Sersic2D(\n", "\n", " amplitude=0.1, # Intensity at r_eff\n", " r_eff=10, # Effective or half-lilght radius\n", " n=1.7384901, # Sersic index\n", " x_0=99.97722657736085, # center of model in the x direction\n", " y_0=99.12324178530918, # center of model in the y direction\n", " ellip=0.1, # Ellipticity\n", " theta=0.0, # Rotation angle in radians, counterclockwise from the positive x-axis.\n", "\n", " bounds=pf.get_default_sersic_bounds(), # Parameter bounds\n", ") \n", "\n", "# Football shaped galaxy \n", "galaxy_model_2 = models.Sersic2D(\n", "\n", " amplitude=0.1, # Intensity at r_eff\n", " r_eff=10, # Effective or half-lilght radius\n", " n=1, # Sersic index\n", " x_0=138.56315299695075, # center of model in the x direction\n", " y_0=89.27757468116197, # center of model in the y direction\n", " ellip=0.7, # Ellipticity\n", " theta=0.7, # Rotation angle in radians, counterclockwise from the positive x-axis.\n", "\n", " bounds=pf.get_default_sersic_bounds(), # Parameter bounds\n", ")\n", "\n", "# Large galaxy near the bottom corner \n", "galaxy_model_3 = models.Sersic2D(\n", "\n", " amplitude=0.1, # Intensity at r_eff\n", " r_eff=10, # Effective or half-lilght radius\n", " n=1, # Sersic index\n", " x_0=178.72302596615611, # center of model in the x direction\n", " y_0=63.506754312433046\t, # center of model in the y direction\n", " ellip=0.2, # Ellipticity\n", " theta=0.0, # Rotation angle in radians, counterclockwise from the positive x-axis.\n", "\n", " bounds=pf.get_default_sersic_bounds(), # Parameter bounds\n", ") " ] }, { "cell_type": "markdown", "id": "b211ac6c", "metadata": {}, "source": [ "Make compound PSF model as follows:" ] }, { "cell_type": "code", "execution_count": null, "id": "780fef6e", "metadata": {}, "outputs": [], "source": [ "compound_galaxies_model = galaxy_model_1 + galaxy_model_2 + galaxy_model_3\n", "\n", "all_galaxies_psf_model = pf.PSFConvolvedModel2D(compound_galaxies_model, psf=PSF)" ] }, { "cell_type": "markdown", "id": "4c368d95", "metadata": {}, "source": [ "Fit the model " ] }, { "cell_type": "code", "execution_count": null, "id": "5fa15fd3", "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "all_galaxies_fitted_model, fitter = pf.fit_model(\n", " image=image.data, \n", " model=all_galaxies_psf_model,\n", " weights=1/rms, # optional \n", " maxiter=10000,\n", " epsilon=1.4901161193847656e-08,\n", " acc=1e-09,\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "bb2760ad", "metadata": {}, "outputs": [], "source": [ "pf.plot_fit(all_galaxies_fitted_model, image, vmax=vmax, vmin=vmin, figsize=[6*3, 6])\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "10207acc", "metadata": {}, "source": [ "Looks like the bottom corner galaxy is a spiral, let us add another component for the spiral and use `LevMarLSQFitter` to fit again:" ] }, { "cell_type": "code", "execution_count": null, "id": "560572e8", "metadata": {}, "outputs": [], "source": [ "%%time\n", "from astropy.modeling import fitting\n", "\n", "# Redefine model with an extra component for galaxy 3 \n", "compound_galaxies_model = galaxy_model_1 + galaxy_model_2 + galaxy_model_3 + galaxy_model_3\n", "\n", "# PSF model\n", "all_galaxies_psf_model = pf.PSFConvolvedModel2D(compound_galaxies_model, psf=PSF)\n", "\n", "# Fit the model\n", "all_galaxies_fitted_model, fitter = pf.fit_model(\n", " image=image.data, \n", " model=all_galaxies_psf_model,\n", " weights=1/rms, # optional\n", " maxiter=100000,\n", " fitter=fitting.LevMarLSQFitter,\n", " epsilon=1.4901161193847656e-08,\n", " acc=1e-09,\n", ")\n", "\n", "# Plot the fit\n", "pf.plot_fit(all_galaxies_fitted_model, image, vmax=vmax, vmin=vmin, figsize=[6*3, 6])\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "ed638cd8", "metadata": {}, "source": [ "## Fitting Image Backgrounds \n", "\n", "The `fit_background` function can be used to fit the background pixels using a 2D plane. It will sigma clip the pixels (sigma value provided by the user) and fit a 2D plane to the clipped image. Users can also provide their own 2D models." ] }, { "cell_type": "code", "execution_count": null, "id": "0a7f5fb5", "metadata": {}, "outputs": [], "source": [ "bg_model, fitter = pf.fit_background(image, sigma=3.0)\n", "bg_image = pf.model_to_image(bg_model, size=(image.shape[1], image.shape[0]))" ] }, { "cell_type": "markdown", "id": "6cef5365", "metadata": {}, "source": [ "Plot backround" ] }, { "cell_type": "code", "execution_count": null, "id": "d664f153", "metadata": {}, "outputs": [], "source": [ "fig, axs = plt.subplots(1,2, figsize=[6*2, 6])\n", "\n", "plt.sca(axs[0])\n", "plt.imshow(bg_image)\n", "plt.title(\"Background Image\")\n", "plt.xlabel(\"Pixels\")\n", "plt.ylabel(\"Pixels\")\n", "\n", "plt.sca(axs[1])\n", "plt.imshow(image.data - bg_image, vmin=vmin, vmax=vmax)\n", "plt.title(\"Background subtracted image\")\n", "plt.xlabel(\"Pixels\")\n", "plt.ylabel(\"Pixels\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "a9df1d7b", "metadata": {}, "source": [ "## Fitting a PSF With a Moffat Model\n", "\n", "In this example we fit the PSF itself, from [the section above](#Making-a-PSF-Convolved-Model), using an `astropy` 2D Moffat model (PSF convolution is not needed for such fits). We then use the model PSF to fit the cutout image from the [Making Cutouts section](#Making-Cutouts). We start by initializing a `Moffat2D` model:" ] }, { "cell_type": "code", "execution_count": null, "id": "fbe03244", "metadata": {}, "outputs": [], "source": [ "moffat_model = models.Moffat2D(amplitude=1, x_0=25, y_0=25, gamma=1, alpha=1)" ] }, { "cell_type": "markdown", "id": "693f7dae", "metadata": {}, "source": [ "Fit the model using ``fit_model``:" ] }, { "cell_type": "code", "execution_count": null, "id": "0b8fa2d5", "metadata": {}, "outputs": [], "source": [ "%%time\n", "\n", "fitted_moffat_model, fitter = pf.fit_model(\n", " image=PSF, \n", " model=moffat_model,\n", " weights=None, # optional \n", " maxiter=10000,\n", " epsilon=1.4901161193847656e-08,\n", " acc=1e-09,\n", ")" ] }, { "cell_type": "markdown", "id": "5aaea2cc", "metadata": {}, "source": [ "Plot the fit and print out fitted parameters:" ] }, { "cell_type": "code", "execution_count": null, "id": "db230d62", "metadata": {}, "outputs": [], "source": [ "psf_vmax = PSF.std()/10\n", "psf_vmin = -psf_vmax\n", "pf.plot_fit(fitted_moffat_model, PSF, vmin=psf_vmin, vmax=psf_vmax, figsize=[6*3, 6])\n", "plt.show()\n", "\n", "print(\"Fitted Moffat Params:\")\n", "pf.print_model_params(fitted_moffat_model)" ] }, { "cell_type": "markdown", "id": "e002d052", "metadata": {}, "source": [ "Use the fitted Moffat model as a PSF and fit a galaxy:" ] }, { "cell_type": "code", "execution_count": null, "id": "3db0accd", "metadata": {}, "outputs": [], "source": [ "# Make Moffat PSF \n", "moffat_model_psf = pf.model_to_image(fitted_moffat_model, size=(51, 51))\n", "\n", "# Make a PSFConvolvedModel2D model with Moffat PSF \n", "moffat_psf_sersic_model = pf.PSFConvolvedModel2D(\n", " sersic_model, \n", " psf=moffat_model_psf, # Moffat PSF \n", " oversample=4\n", ")\n", "\n", "# Fit the galaxy cutout image\n", "fitted_moffat_psf_sersic_model, fitter = pf.fit_model(\n", " image=cutout_image.data, \n", " model=moffat_psf_sersic_model,\n", " weights=fitting_weights, # optional \n", " maxiter=10000,\n", " epsilon=1.4901161193847656e-08,\n", " acc=1e-09,\n", ")\n", "\n", "pf.plot_fit(fitted_moffat_psf_sersic_model, cutout_image.data, vmin=vmin, vmax=vmax, figsize=[6*3, 6])\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "d93fc8d5", "metadata": {}, "source": [ "## Making a Photutils Source Catalog\n", "\n", "To make a Photutils source catalog, which can also be converted into a table, use the `make_catalog` wrapper as follows:" ] }, { "cell_type": "code", "execution_count": null, "id": "4e459af0", "metadata": {}, "outputs": [], "source": [ "from astropy.stats import sigma_clipped_stats\n", "\n", "# Sigma clipped stats\n", "image_mean, image_median, image_stddev = sigma_clipped_stats(image.data, sigma=3)\n", "\n", "cat, segm, segm_deblend = pf.make_catalog(\n", " image=image.data, # Input image\n", " threshold=image_stddev*3, # Detection threshold\n", " deblend=True, # Deblend sources?\n", " npixels=4**2, # Minimum number of pixels that make up a source\n", " plot=True, vmax=vmax, vmin=vmin, # Plotting params\n", " figsize=(12, 6)\n", ")" ] }, { "cell_type": "markdown", "id": "62803414", "metadata": {}, "source": [ "Display sources catalog as table:" ] }, { "cell_type": "code", "execution_count": null, "id": "9d638ef2", "metadata": {}, "outputs": [], "source": [ "cat.to_table()" ] }, { "cell_type": "markdown", "id": "3f83da65", "metadata": {}, "source": [ "## Curve of Growth and Petrosian Radii\n", "\n", "This step in detailed in the [Photometry](./photometry.ipynb#Photometry) and [Petrosian](./petrosian.ipynb#Petrosian) sections but we will do a simple example here. The first step is to pick a source from the catalog we made in the last step:" ] }, { "cell_type": "code", "execution_count": null, "id": "25f08e87", "metadata": {}, "outputs": [], "source": [ "source = cat[6]\n", "\n", "# Photutils cutout of the source\n", "# Not to be confused with the cutout we made\n", "plt.imshow(source.data, vmin=vmin, vmax=vmax)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "b56fb3d1", "metadata": {}, "source": [ "Now we use PetroFit `source_photometry` to compute the curve of growth" ] }, { "cell_type": "code", "execution_count": null, "id": "bc926566", "metadata": {}, "outputs": [], "source": [ "r_list = pf.make_radius_list(\n", " max_pix=50, # Max pixel to go up to\n", " n=50 # the number of radii to produce\n", ")\n", "\n", "# Photomerty\n", "flux_arr, area_arr, error_arr = pf.source_photometry(\n", "\n", " # Inputs\n", " source, # Source (`photutils.segmentation.catalog.SourceCatalog`)\n", " image.data, # Image as 2D array\n", " segm_deblend, # Deblended segmentation map of image\n", " r_list, # list of aperture radii\n", "\n", " # Options\n", " cutout_size=max(r_list)*2, # Cutout out size, set to double the max radius\n", " bg_sub=True, # Subtract background\n", " sigma=3, sigma_type='clip', # Fit a 2D plane to pixels within 3 sigma of the mean\n", " plot=True, vmax=vmax, vmin=vmin, # Show plot with max and min defined above\n", ")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "24552d18", "metadata": {}, "source": [ "Now we have a radius list (`r_list`) and a corresponding enclosed flux list (`flux_arr`), we can plot the curve of growth and initialize a `Petrosian` profile" ] }, { "cell_type": "code", "execution_count": null, "id": "c83c5aa8", "metadata": {}, "outputs": [], "source": [ "plt.plot(r_list, flux_arr, linewidth=3, marker='o')\n", "\n", "pf.mpl_tick_frame()\n", "\n", "plt.title(\"Curve of Growth\")\n", "plt.xlabel(\"Radius R [Pix]\")\n", "plt.ylabel(\"$L(