Add Source Catalog Step

CHANGES.rst

@@ -92,6 +92,11 @@ flux
 
 - Create FluxStep to apply the flux correction to Level 2 data. [#1120]
 
+source_catalog
+--------------
+
+- Added Source Catalog Step. [#1102]
+
 
 0.14.0 (2024-02-12)
 ===================

docs/roman/package_index.rst

@@ -16,6 +16,7 @@
    refpix/index.rst
    resample/index.rst
    skymatch/index.rst
+   source_catalog/index.rst
    source_detection/index.rst
    stpipe/index.rst
    tweakreg/index.rst

docs/roman/pipeline_steps.rst

@@ -19,3 +19,4 @@ Pipeline Steps
    skymatch/index.rst
    outlier_detection/index.rst
    resample/index.rst
+   source_catalog/index.rst

docs/roman/source_catalog/arguments.rst

@@ -0,0 +1,43 @@
+Arguments
+=========
+
+The ``source_catalog`` step uses the following user-settable arguments:
+
+* ``--bkg_boxsize``: An integer value giving the background mesh box
+  size in pixels
+
+* ``--kernel_fwhm``: A floating-point value giving the Gaussian kernel
+  FWHM in pixels
+
+* ``--snr_threshold``: A floating-point value that sets the SNR
+  threshold (above background) for source detection
+
+* ``--npixels``: An integer value that sets the minimum number of
+  pixels in a source
+
+* ``--deblend``: A boolean indicating whether to deblend sources
+
+
+* ``--aperture_ee1``: An integer value of the smallest aperture
+  encircled energy value
+
+* ``--aperture_ee2``: An integer value of the middle aperture encircled
+  energy value
+
+* ``--aperture_ee3``: An integer value of the largest aperture encircled
+  energy value
+
+* ``--ci1_star_threshold``: A floating-point value of the threshold
+  compared to the concentration index calculated from aperture_ee1
+  and aperture_ee2 that is used to determine whether a source is a
+  star. Sources must meet the criteria of both ci1_star_threshold and
+  ci2_star_threshold to be considered a star.
+
+* ``--ci2_star_threshold``: A floating-point value of the threshold
+  compared to the concentration index calculated from aperture_ee2
+  and aperture_ee3 that is used to determine whether a source is a
+  star. Sources must meet the criteria of both ci1_star_threshold and
+  ci2_star_threshold to be considered a star.
+
+* ``--suffix``: A string value giving the file name suffix to use for
+  the output catalog file [default='cat']

docs/roman/source_catalog/index.rst

@@ -0,0 +1,13 @@
+.. _source_catalog_step:
+
+==============
+Source Catalog
+==============
+
+.. toctree::
+   :maxdepth: 2
+
+   main.rst
+   arguments.rst
+
+.. automodapi:: romancal.source_catalog

docs/roman/source_catalog/main.rst

@@ -0,0 +1,184 @@
+Description
+===========
+
+:Class: `romancal.source_catalog.SourceCatalogStep`
+:Alias: source_catalog
+
+This step creates a catalog of source photometry and morphologies.
+Both aperture and isophotal (segment-based) photometry are calculated.
+Source morphologies are based on 2D image moments within the source
+segment.
+
+
+Source Detection
+----------------
+Sources are detected using `image segmentation
+<https://en.wikipedia.org/wiki/Image_segmentation>`_, which is a
+process of assigning a label to every pixel in an image such that
+pixels with the same label are part of the same source.  The
+segmentation procedure used is from `Photutils source extraction
+<https://photutils.readthedocs.io/en/latest/segmentation.html>`_.
+Detected sources must have a minimum number of connected pixels that
+are each greater than a specified threshold value in an image.  The
+threshold level is usually defined at some multiple of the background
+standard deviation above the background.  The image can also be
+filtered before thresholding to smooth the noise and maximize the
+detectability of objects with a shape similar to the filter kernel.
+
+Source Deblending
+-----------------
+Overlapping sources are detected as single sources.  Separating those
+sources requires a deblending procedure, such as a multi-thresholding
+technique used by `SExtractor
+<https://www.astromatic.net/software/sextractor>`_.  Here we use the
+`Photutils deblender
+<https://photutils.readthedocs.io/en/latest/segmentation.html#source-deblending>`_,
+which is an algorithm that deblends sources using a combination of
+multi-thresholding and `watershed segmentation
+<https://en.wikipedia.org/wiki/Watershed_(image_processing)>`_.  In
+order to deblend sources, they must be separated enough such that
+there is a saddle between them.
+
+Source Photometry and Properties
+--------------------------------
+After detecting sources using image segmentation, we can measure their
+photometry, centroids, and morphological properties.  The aperture
+photometry is measured in three apertures, based on the input
+encircled energy values.  The total aperture-corrected flux and
+magnitudes are also calculated, based on the largest aperture.
+
+The isophotal photometry is based on `photutils segmentation
+<https://photutils.readthedocs.org/en/latest/segmentation.html>`_.
+The properties that are currently calculated for each source include
+source centroids (both in pixel and sky coordinates), isophotal fluxes
+(and errors), isophotal area,
+semimajor and semiminor axis lengths, orientation of the major axis,
+and sky coordinates at corners of the minimal bounding box enclosing
+the source.
+
+Photometric errors are calculated from the resampled total-error
+array contained in the ``ERR`` (``model.err``) array. Note that this
+total-error array includes source Poisson noise.
+
+Output Products
+---------------
+
+Source Catalog Table
+^^^^^^^^^^^^^^^^^^^^
+The output source catalog table is saved in `ECSV format
+<https://docs.astropy.org/en/stable/io/ascii/ecsv.html>`_.
+
+The table contains a row for each source, with the following default
+columns (assuming the default encircled energies of 30, 50, and 70):
+
++------------------------+----------------------------------------------------+
+| Column                 | Description                                        |
++========================+====================================================+
+| label                  | Unique source identification label number          |
++------------------------+----------------------------------------------------+
+| xcentroid              | X pixel value of the source centroid (0 indexed)   |
++------------------------+----------------------------------------------------+
+| ycentroid              | Y pixel value of the source centroid (0 indexed)   |
++------------------------+----------------------------------------------------+
+| sky_centroid           | Sky coordinate of the source centroid              |
++------------------------+----------------------------------------------------+
+| aper_bkg_flux          | The local background value calculated as the       |
+|                        | sigma-clipped median value in the background       |
+|                        | annulus aperture                                   |
++------------------------+----------------------------------------------------+
+| aper_bkg_flux_err      | The standard error of the sigma-clipped median     |
+|                        | background value                                   |
++------------------------+----------------------------------------------------+
+| aper30_flux            | Flux within the 30% encircled energy circular      |
+|                        | aperture                                           |
++------------------------+----------------------------------------------------+
+| aper30_flux_err        | Flux error within the 30% encircled energy         |
+|                        | circular aperture                                  |
++------------------------+----------------------------------------------------+
+| aper50_flux            | Flux within the 50% encircled energy circular      |
+|                        | aperture                                           |
++------------------------+----------------------------------------------------+
+| aper50_flux_err        | Flux error within the 50% encircled energy         |
+|                        | circular aperture                                  |
++------------------------+----------------------------------------------------+
+| aper70_flux            | Flux within the 70% encircled energy circular      |
+|                        | aperture                                           |
++------------------------+----------------------------------------------------+
+| aper70_flux_err        | Flux error within the 70% encircled energy         |
+|                        | circular aperture                                  |
++------------------------+----------------------------------------------------+
+| aper_total_flux        | Total aperture-corrected flux based on the 70%     |
+|                        | encircled energy circular aperture; should be used |
+|                        | only for unresolved sources                        |
++------------------------+----------------------------------------------------+
+| aper_total_flux_err    | Total aperture-corrected flux error based on the   |
+|                        | 70% encircled energy circular aperture; should be  |
+|                        | used only for unresolved sources                   |
++------------------------+----------------------------------------------------+
+| CI_50_30               | Concentration index calculated as (aper50_flux /   |
+|                        | aper30_flux)                                       |
++------------------------+----------------------------------------------------+
+| CI_70_50               | Concentration index calculated as (aper70_flux /   |
+|                        | aper50_flux)                                       |
++------------------------+----------------------------------------------------+
+| CI_70_30               | Concentration index calculated as (aper70_flux /   |
+|                        | aper30_flux)                                       |
++------------------------+----------------------------------------------------+
+| is_extended            | Flag indicating whether the source is extended     |
++------------------------+----------------------------------------------------+
+| sharpness              | The DAOFind source sharpness statistic             |
++------------------------+----------------------------------------------------+
+| roundness              | The DAOFind source roundness statistic             |
++------------------------+----------------------------------------------------+
+| nn_label               | The label number of the nearest neighbor           |
++------------------------+----------------------------------------------------+
+| nn_dist                | The distance in pixels to the nearest neighbor     |
++------------------------+----------------------------------------------------+
+| isophotal_flux         | Isophotal flux                                     |
++------------------------+----------------------------------------------------+
+| isophotal_flux_err     | Isophotal flux error                               |
++------------------------+----------------------------------------------------+
+| isophotal_area         | Isophotal area                                     |
++------------------------+----------------------------------------------------+
+| semimajor_sigma        | 1-sigma standard deviation along the semimajor     |
+|                        | axis of the 2D Gaussian function that has the same |
+|                        | second-order central moments as the source         |
++------------------------+----------------------------------------------------+
+| semiminor_sigma        | 1-sigma standard deviation along the semiminor     |
+|                        | axis of the 2D Gaussian function that has the same |
+|                        | second-order central moments as the source         |
++------------------------+----------------------------------------------------+
+| ellipticity            | 1 minus the ratio of the 1-sigma lengths of the    |
+|                        | semimajor and semiminor axes                       |
++------------------------+----------------------------------------------------+
+| orientation            | The angle (degrees) between the positive X axis    |
+|                        | and the major axis (increases counter-clockwise)   |
++------------------------+----------------------------------------------------+
+| sky_orientation        | The position angle (degrees) from North of the     |
+|                        | major axis                                         |
++------------------------+----------------------------------------------------+
+| sky_bbox_ll            | Sky coordinate of the lower-left vertex of the     |
+|                        | minimal bounding box of the source                 |
++------------------------+----------------------------------------------------+
+| sky_bbox_ul            | Sky coordinate of the upper-left vertex of the     |
+|                        | minimal bounding box of the source                 |
++------------------------+----------------------------------------------------+
+| sky_bbox_lr            | Sky coordinate of the lower-right vertex of the    |
+|                        | minimal bounding box of the source                 |
++------------------------+----------------------------------------------------+
+| sky_bbox_ur            | Sky coordinate of the upper-right vertex of the    |
+|                        | minimal bounding box of the source                 |
++------------------------+----------------------------------------------------+
+
+Note that pixel coordinates are 0 indexed, matching the Python 0-based
+indexing. That means pixel coordinate ``0`` is the center of the first
+pixel.
+
+
+Segmentation Map
+^^^^^^^^^^^^^^^^
+
+The segmentation map computed during the source finding process is saved
+to a single 2D image extension in a FITS file. Each image pixel contains an
+integer value corresponding to a source label number in the source catalog
+product. Pixels that don't belong to any source have a value of zero.

romancal/lib/suffix.py

@@ -60,6 +60,9 @@
     "outlier_detection",
     "skymatch",
     "refpix",
+    "cat",
+    "segm",
+    "i2d",
 ]
 
 # Suffixes that are discovered but should not be considered.
@@ -95,6 +98,7 @@
     "skymatchstep",
     "refpixstep",
     "resamplestep",
+    "sourcecatalogstep",
 }
 
 

romancal/regtest/test_catalog.py

@@ -0,0 +1,102 @@
+""" Roman tests for source catalog creation """
+
+import asdf
+import pytest
+from metrics_logger.decorators import metrics_logger
+
+from romancal.source_catalog.source_catalog_step import SourceCatalogStep
+from romancal.stpipe import RomanStep
+
+
+def passfail(bool_expr):
+    """set pass fail"""
+    if bool_expr:
+        return "Pass"
+    return "Fail"
+
+
+def check_catalog_fields(model, log, modeltype):
+    fields = model.dtype.names
+    has_pos = ("ra_centroid" in fields) and ("dec_centroid" in fields)
+    log.info(
+        f"DMS374 MSG: {modeltype} Catalog contains sky coordinates of sources? :"
+        + passfail(has_pos)
+    )
+
+    has_flux = "aper_total_flux" in fields
+    log.info(f"DMS375 MSG: {modeltype} Catalog contains fluxes? :" + passfail(has_flux))
+
+    has_type = "is_extended" in fields
+    log.info(
+        f"DMS376 MSG: {modeltype} Catalog contains source classification? :"
+        + passfail(has_type)
+    )
+
+    has_flux_err = "aper_total_flux_err" in fields
+    log.info(
+        f"DMS386 MSG: {modeltype} Catalog contains flux uncertainties? :"
+        + passfail(has_flux_err)
+    )
+
+    has_flags = "flags" in fields
+    log.info(
+        f"DMS387 MSG: {modeltype} Catalog contains DQ flags? :" + passfail(has_flags)
+    )
+
+    return has_pos & has_flux & has_type & has_flux_err & has_flags
+
+
+@pytest.mark.bigdata
+@pytest.mark.soctests
+@metrics_logger("DMS374", "DMS375", "DMS376", "DMS386", "DMS387")
+def test_catalog_l3(rtdata, ignore_asdf_paths):
+    # DMS374: positions on ICRF
+    # DMS375: fluxes
+    # DMS376: type of source
+    # DMS386: flux uncertainties
+    # DMS387: DQ flags
+    inputfn = "r0099101001001001001_F158_visit_0.900.0.50_178199.5_-0.5_i2d.asdf"
+    outputfn = inputfn.replace("_i2d", "_cat")
+    rtdata.get_data(f"WFI/image/{inputfn}")
+    rtdata.input = inputfn
+    rtdata.get_truth(f"truth/WFI/image/{outputfn}")
+
+    args = [
+        "romancal.step.SourceCatalogStep",
+        rtdata.input,
+    ]
+    RomanStep.from_cmdline(args)
+    catalogfp = asdf.open(outputfn)
+    catalog = catalogfp["roman"]["source_catalog"]
+    step = SourceCatalogStep()
+    assert check_catalog_fields(catalog, step.log, "L3")
+
+    # no compare_asdf on the catalogs
+
+
+@pytest.mark.bigdata
+@pytest.mark.soctests
+@metrics_logger("DMS374", "DMS375", "DMS376", "DMS386", "DMS387")
+def test_catalog_l2(rtdata, ignore_asdf_paths):
+    # DMS374: positions on ICRF
+    # DMS375: fluxes
+    # DMS376: type of source
+    # DMS386: flux uncertainties
+    # DMS387: DQ flags
+    inputfn = "r0000101001001001001_01101_0001_WFI01_cal.asdf"
+    outputfn = inputfn.replace("cal", "cat")
+    rtdata.get_data(f"WFI/image/{inputfn}")
+    rtdata.input = inputfn
+    rtdata.get_truth(f"truth/WFI/image/{outputfn}")
+
+    args = [
+        "romancal.step.SourceCatalogStep",
+        rtdata.input,
+    ]
+    RomanStep.from_cmdline(args)
+    catalogfp = asdf.open(outputfn)
+    catalog = catalogfp["roman"]["source_catalog"]
+    step = SourceCatalogStep()
+
+    assert check_catalog_fields(catalog, step.log, "L2")
+    # no compare_asdf on the catalogs