From a25180ec0927e142f192ecaf49d697fb12cb3463 Mon Sep 17 00:00:00 2001 From: "Brett M. Morris" Date: Mon, 12 Feb 2024 13:17:26 -0500 Subject: [PATCH] add source catalog step --- docs/roman/package_index.rst | 1 + docs/roman/pipeline_steps.rst | 1 + docs/roman/source_catalog/arguments.rst | 43 + docs/roman/source_catalog/index.rst | 13 + docs/roman/source_catalog/main.rst | 184 +++ romancal/source_catalog/__init__.py | 3 + romancal/source_catalog/_wcs_helpers.py | 65 ++ romancal/source_catalog/detection.py | 254 ++++ romancal/source_catalog/reference_data.py | 264 +++++ romancal/source_catalog/source_catalog.py | 1036 +++++++++++++++++ .../source_catalog/source_catalog_step.py | 123 ++ romancal/source_catalog/tests/__init__.py | 0 .../tests/test_source_catalog.py | 159 +++ 13 files changed, 2146 insertions(+) create mode 100644 docs/roman/source_catalog/arguments.rst create mode 100644 docs/roman/source_catalog/index.rst create mode 100644 docs/roman/source_catalog/main.rst create mode 100644 romancal/source_catalog/__init__.py create mode 100644 romancal/source_catalog/_wcs_helpers.py create mode 100644 romancal/source_catalog/detection.py create mode 100644 romancal/source_catalog/reference_data.py create mode 100644 romancal/source_catalog/source_catalog.py create mode 100755 romancal/source_catalog/source_catalog_step.py create mode 100644 romancal/source_catalog/tests/__init__.py create mode 100644 romancal/source_catalog/tests/test_source_catalog.py diff --git a/docs/roman/package_index.rst b/docs/roman/package_index.rst index ddf4b70e5..efa1ae0aa 100644 --- a/docs/roman/package_index.rst +++ b/docs/roman/package_index.rst @@ -15,6 +15,7 @@ refpix/index.rst resample/index.rst skymatch/index.rst + source_catalog/index.rst source_detection/index.rst stpipe/index.rst tweakreg/index.rst diff --git a/docs/roman/pipeline_steps.rst b/docs/roman/pipeline_steps.rst index 8c44e8cba..33e0214f5 100644 --- a/docs/roman/pipeline_steps.rst +++ b/docs/roman/pipeline_steps.rst @@ -12,6 +12,7 @@ Pipeline Steps assign_wcs/index.rst flatfield/index.rst photom/index.rst + source_catalog/index.rst source_detection/index.rst tweakreg/index.rst resample/index.rst diff --git a/docs/roman/source_catalog/arguments.rst b/docs/roman/source_catalog/arguments.rst new file mode 100644 index 000000000..26883c43e --- /dev/null +++ b/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'] diff --git a/docs/roman/source_catalog/index.rst b/docs/roman/source_catalog/index.rst new file mode 100644 index 000000000..6c1221ee8 --- /dev/null +++ b/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 diff --git a/docs/roman/source_catalog/main.rst b/docs/roman/source_catalog/main.rst new file mode 100644 index 000000000..fea583d8e --- /dev/null +++ b/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 +`_, 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 +`_. +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 +`_. Here we use the +`Photutils deblender +`_, +which is an algorithm that deblends sources using a combination of +multi-thresholding and `watershed segmentation +`_. 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 +`_. +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 +`_. + +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. diff --git a/romancal/source_catalog/__init__.py b/romancal/source_catalog/__init__.py new file mode 100644 index 000000000..faf55c960 --- /dev/null +++ b/romancal/source_catalog/__init__.py @@ -0,0 +1,3 @@ +from .source_catalog_step import SourceCatalogStep + +__all__ = ["SourceCatalogStep"] diff --git a/romancal/source_catalog/_wcs_helpers.py b/romancal/source_catalog/_wcs_helpers.py new file mode 100644 index 000000000..59b5e1551 --- /dev/null +++ b/romancal/source_catalog/_wcs_helpers.py @@ -0,0 +1,65 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +# (taken from photutils: should probably migrate into astropy.wcs) +""" +This module provides WCS helper tools. +""" + +import astropy.units as u +import numpy as np + + +def pixel_scale_angle_at_skycoord(skycoord, wcs, offset=1 * u.arcsec): + """ + Calculate the pixel coordinate and scale and WCS rotation angle at + the position of a SkyCoord coordinate. + + Parameters + ---------- + skycoord : `~astropy.coordinates.SkyCoord` + The SkyCoord coordinate. + + wcs : WCS object + A world coordinate system (WCS) transformation that + supports the `astropy shared interface for WCS + `_ (e.g., + `astropy.wcs.WCS`, `gwcs.wcs.WCS`). + + offset : `~astropy.units.Quantity` + A small angular offset to use to compute the pixel scale and + position angle. + + Returns + ------- + xypos : tuple of float + The (x, y) pixel coordinate. + + scale : `~astropy.units.Quantity` + The pixel scale in arcsec/pixel. + + angle : `~astropy.units.Quantity` + The angle (in degrees) measured counterclockwise from the + positive x axis to the "North" axis of the celestial coordinate + system. + + Notes + ----- + If distortions are present in the image, the x and y pixel scales + likely differ. This function computes a single pixel scale along + the North/South axis. + """ + # Convert to pixel coordinates + xpos, ypos = wcs.world_to_pixel(skycoord) + + # We take a point directly North (i.e., latitude offset) the + # input sky coordinate and convert it to pixel coordinates, + # then we use the pixel deltas between the input and offset sky + # coordinate to calculate the pixel scale and angle. + skycoord_offset = skycoord.directional_offset_by(0.0, offset) + x_offset, y_offset = wcs.world_to_pixel(skycoord_offset) + + dx = x_offset - xpos + dy = y_offset - ypos + scale = offset.to(u.arcsec) / (np.hypot(dx, dy) * u.pixel) + angle = (np.arctan2(dy, dx) * u.radian).to(u.deg) + + return (xpos, ypos), scale, angle diff --git a/romancal/source_catalog/detection.py b/romancal/source_catalog/detection.py new file mode 100644 index 000000000..795ba5756 --- /dev/null +++ b/romancal/source_catalog/detection.py @@ -0,0 +1,254 @@ +""" +Module to detect sources using image segmentation. +""" + +import logging +import warnings + +import numpy as np +from astropy.convolution import Gaussian2DKernel, convolve +from astropy.stats import SigmaClip, gaussian_fwhm_to_sigma +from astropy.utils import lazyproperty +from astropy.utils.exceptions import AstropyUserWarning +from photutils.background import Background2D, MedianBackground +from photutils.segmentation import deblend_sources, detect_sources +from photutils.utils.exceptions import NoDetectionsWarning + +log = logging.getLogger(__name__) +log.setLevel(logging.DEBUG) + + +class RomanBackground: + """ + Class to estimate a 2D background and background RMS noise in an + image. + + Parameters + ---------- + data : 2D `~numpy.ndarray` + The input 2D array. + + box_size : int or array_like (int) + The box size along each axis. If ``box_size`` is a scalar then + a square box of size ``box_size`` will be used. If ``box_size`` + has two elements, they should be in ``(ny, nx)`` order. + + coverage_mask : array_like (bool), optional + A boolean mask, with the same shape as ``data``, where a `True` + value indicates the corresponding element of ``data`` is masked. + Masked data are excluded from calculations. ``coverage_mask`` + should be `True` where there is no coverage (i.e., no data) for + a given pixel (e.g., blank areas in a mosaic image). It should + not be used for bad pixels. + + Attributes + ---------- + background : 2D `~numpy.ndimage` + The estimated 2D background image. + + background_rms : 2D `~numpy.ndimage` + The estimated 2D background RMS image. + """ + + def __init__(self, data, box_size=100, coverage_mask=None): + self.data = data + self.box_size = np.asarray(box_size).astype(int) # must be integer + self.coverage_mask = coverage_mask + + @lazyproperty + def _background2d(self): + """ + Estimate the 2D background and background RMS noise in an image. + + Returns + ------- + background : `photutils.background.Background2D` + A Background2D object containing the 2D background and + background RMS noise estimates. + """ + sigma_clip = SigmaClip(sigma=3.0) + bkg_estimator = MedianBackground() + filter_size = (3, 3) + + try: + bkg = Background2D( + self.data, + self.box_size, + filter_size=filter_size, + coverage_mask=self.coverage_mask, + sigma_clip=sigma_clip, + bkg_estimator=bkg_estimator, + ) + except ValueError: + # use the entire unmasked array + bkg = Background2D( + self.data, + self.data.shape, + filter_size=filter_size, + coverage_mask=self.coverage_mask, + sigma_clip=sigma_clip, + bkg_estimator=bkg_estimator, + exclude_percentile=100.0, + ) + log.info( + "Background could not be estimated in meshes. " + "Using the entire unmasked array for background " + f"estimation: bkg_boxsize={self.data.shape}." + ) + + return bkg + + @lazyproperty + def background(self): + """ + The 2D background image. + """ + return self._background2d.background + + @lazyproperty + def background_rms(self): + """ + The 2D background RMS image. + """ + return self._background2d.background_rms + + +def make_kernel(kernel_fwhm): + """ + Make a 2D Gaussian smoothing kernel that is used to filter the image + before thresholding. + + Filtering the image will smooth the noise and maximize detectability + of objects with a shape similar to the kernel. + + The kernel must have odd sizes in both X and Y, be centered in the + central pixel, and normalized to sum to 1. + + Parameters + ---------- + kernel_fwhm : float + The full-width at half-maximum (FWHM) of the 2D Gaussian kernel. + + Returns + ------- + kernel : `astropy.convolution.Kernel2D` + The output smoothing kernel, normalized such that it sums to 1. + """ + sigma = kernel_fwhm * gaussian_fwhm_to_sigma + kernel = Gaussian2DKernel(sigma) + kernel.normalize(mode="integral") + return kernel + + +def convolve_data(data, kernel_fwhm, mask=None): + """ + Convolve the data with a Gaussian2D kernel. + + Parameters + ---------- + data : `~numpy.ndarray` + The 2D array to convolve. + + kernel_fwhm : float + The full-width at half-maximum (FWHM) of the 2D Gaussian kernel. + + mask : array_like, bool, optional + A boolean mask with the same shape as ``data``, where a `True` + value indicates the corresponding element of ``data`` is masked. + """ + kernel = make_kernel(kernel_fwhm) + + with warnings.catch_warnings(): + warnings.simplefilter("ignore", AstropyUserWarning) + return convolve(data, kernel, mask=mask, normalize_kernel=True) + + +class RomanSourceFinder: + """ + Class to detect sources, including deblending, using image + segmentation. + + Parameters + ---------- + threshold : float + The data value to be used as the per-pixel detection threshold. + + npixels : int + The number of connected pixels, each greater than the threshold, + that an object must have to be detected. ``npixels`` must be a + positive integer. + + deblend : bool, optional + Whether to deblend overlapping sources. Source deblending + requires scikit-image. + """ + + def __init__(self, threshold, npixels, deblend=False): + self.threshold = threshold + self.npixels = npixels + self.deblend = deblend + self.connectivity = 8 + self.nlevels = 32 + self.contrast = 0.001 + self.mode = "exponential" + + def __call__(self, convolved_data, mask=None): + """ + Parameters + ---------- + convolved_data : 2D `numpy.ndarray` + The 2D convolved array from which to detect sources. + + mask : array_like, bool, optional + A boolean mask with the same shape as ``convolved_data``, + where a `True` value indicates the corresponding element + of ``convolved_data`` is masked. Masked pixels will not be + included in any source. + + Returns + ------- + segment_image : `~photutils.segmentation.SegmentationImage` or `None` + A 2D segmentation image, with the same shape as the input data, + where sources are marked by different positive integer values. A + value of zero is reserved for the background. If no sources are + found then `None` is returned. + """ + if mask is not None: + if mask.all(): + log.error( + "There are no valid pixels in the image to detect " "sources." + ) + return None + + with warnings.catch_warnings(): + # suppress NoDetectionsWarning from photutils + warnings.filterwarnings("ignore", category=NoDetectionsWarning) + + segment_img = detect_sources( + convolved_data, + self.threshold, + self.npixels, + mask=mask, + connectivity=self.connectivity, + ) + if segment_img is None: + log.warning( + "No sources were found. Source catalog will not " "be created." + ) + return None + + # source deblending requires scikit-image + if self.deblend: + segment_img = deblend_sources( + convolved_data, + segment_img, + npixels=self.npixels, + nlevels=self.nlevels, + contrast=self.contrast, + mode=self.mode, + connectivity=self.connectivity, + relabel=True, + ) + + log.info(f"Detected {segment_img.nlabels} sources") + return segment_img diff --git a/romancal/source_catalog/reference_data.py b/romancal/source_catalog/reference_data.py new file mode 100644 index 000000000..e42de21c1 --- /dev/null +++ b/romancal/source_catalog/reference_data.py @@ -0,0 +1,264 @@ +""" +Module for parsing APCORR and ABVEGAOFFSET reference file data. +""" + +import logging + +import numpy as np +from astropy.utils import lazyproperty +from roman_datamodels import datamodels +from roman_datamodels.datamodels import ImageModel + +log = logging.getLogger(__name__) +log.setLevel(logging.DEBUG) + + +class ReferenceData: + """ + Class for APCORR and ABVEGAOFFSET reference file data needed by + `SourceCatalogStep`. + + Parameters + ---------- + model : `ImageModel` + An `ImageModel` of drizzled image. + + reffile_paths : list of 2 str + The full path filename of the APCORR and ABVEGAOFFSET reference + files. + + aperture_ee : tuple of 3 int + The aperture encircled energies to be used for aperture + photometry. The values must be 3 strictly-increasing integers. + Valid values are defined in the APCORR reference files (20, 30, + 40, 50, 60, 70, or 80). + + Attributes + ---------- + aperture_params : `dict` + A dictionary containing the aperture parameters (radii, aperture + corrections, and background annulus inner and outer radii). + + abvega_offset : float + Offset to convert from AB to Vega magnitudes. The value + represents m_AB - m_Vega. + """ + + def __init__(self, model, reffile_paths, aperture_ee): + if not isinstance(model, ImageModel): + raise ValueError("The input model must be a ImageModel.") + self.model = model + + self.aperture_ee = self._validate_aperture_ee(aperture_ee) + self.apcorr_filename = None + self.abvegaoffset_filename = None + # self.apcorr_filename = reffile_paths[0] + # self.abvegaoffset_filename = reffile_paths[1] + + self.instrument = self.model.meta.instrument.name + self.detector = self.model.meta.instrument.detector + self.filtername = self.model.meta.instrument.filter + # self.pupil = model.meta.instrument.pupil + # self.subarray = self.model.meta.subarray.name + + log.info(f"Instrument: {self.instrument}") + if self.detector is not None: + log.info(f"Detector: {self.detector}") + if self.filtername is not None: + log.info(f"Filter: {self.filtername}") + # if self.pupil is not None: + # log.info(f'Pupil: {self.pupil}') + # if self.subarray is not None: + # log.info(f'Subarray: {self.subarray}') + + @staticmethod + def _validate_aperture_ee(aperture_ee): + """ + Validate the input ``aperture_ee``. + """ + aperture_ee = np.array(aperture_ee).astype(int) + if not np.all(aperture_ee[1:] > aperture_ee[:-1]): + raise ValueError("aperture_ee values must be strictly " "increasing") + if len(aperture_ee) != 3: + raise ValueError("aperture_ee must contain only 3 values") + if np.any(np.logical_or(aperture_ee <= 0, aperture_ee >= 100)): + raise ValueError("aperture_ee values must be between 0 and 100") + return aperture_ee + + @lazyproperty + def _aperture_ee_table(self): + """ + Get the encircled energy table for the given instrument + configuration. + """ + if self.instrument in ("NIRCAM", "NIRISS"): + selector = {"filter": self.filtername, "pupil": self.pupil} + elif self.instrument == "MIRI": + selector = {"filter": self.filtername, "subarray": self.subarray} + elif self.instrument == "FGS": + selector = None + else: + raise RuntimeError(f"{self.instrument} is not a valid instrument") + + apcorr_model = datamodels.open(self.apcorr_filename) + apcorr = apcorr_model.apcorr_table + if selector is None: # FGS + ee_table = apcorr + else: + mask_idx = [apcorr[key] == value for key, value in selector.items()] + ee_table = apcorr[np.logical_and.reduce(mask_idx)] + + if len(ee_table) == 0: + raise RuntimeError( + "APCORR reference file data is missing for " f"{selector}." + ) + + return ee_table + + def _get_ee_table_row(self, aperture_ee): + """ + Get the encircled energy row for the input ``aperture_ee``. + """ + ee_percent = np.round(self._aperture_ee_table["eefraction"] * 100) + row_mask = ee_percent == aperture_ee + ee_row = self._aperture_ee_table[row_mask] + if len(ee_row) == 0: + raise RuntimeError( + "Aperture encircled energy value of " + f"{aperture_ee} appears to be invalid. No " + "matching row was found in the APCORR " + "reference file {self.apcorr_filename}" + ) + if len(ee_row) > 1: + raise RuntimeError( + "More than one matching row was found in " + "the APCORR reference file " + f"{self.apcorr_filename}" + ) + return ee_row + + @lazyproperty + def aperture_params(self): + """ + A dictionary containing the aperture parameters (radii, aperture + corrections, and background annulus inner and outer radii). + """ + if self.apcorr_filename is None: + log.warning( + "APCorrModel reference file was not input. Using " + "fallback aperture sizes without any aperture " + "corrections." + ) + + params = { + "aperture_radii": np.array((1.0, 2.0, 3.0)), + "aperture_corrections": np.array((1.0, 1.0, 1.0)), + "aperture_ee": np.array((1, 2, 3)), + "bkg_aperture_inner_radius": 5.0, + "bkg_aperture_outer_radius": 10.0, + } + return params + + params = {} + radii = [] + apcorrs = [] + skyins = [] + skyouts = [] + for aper_ee in self.aperture_ee: + row = self._get_ee_table_row(aper_ee) + radii.append(row["radius"][0]) + apcorrs.append(row["apcorr"][0]) + skyins.append(row["skyin"][0]) + skyouts.append(row["skyout"][0]) + + if self.model.meta.resample.pixel_scale_ratio is not None: + # pixel_scale_ratio is the ratio of the resampled to the native + # pixel scale (values < 1 have smaller resampled pixels) + pixel_scale_ratio = self.model.meta.resample.pixel_scale_ratio + else: + log.warning( + "model.meta.resample.pixel_scale_ratio was not " + "found. Assuming the native detector pixel scale " + "(i.e., pixel_scale_ratio = 1)" + ) + pixel_scale_ratio = 1.0 + + params["aperture_ee"] = self.aperture_ee + params["aperture_radii"] = np.array(radii) / pixel_scale_ratio + params["aperture_corrections"] = np.array(apcorrs) + + skyins = np.unique(skyins) + skyouts = np.unique(skyouts) + if len(skyins) != 1 or len(skyouts) != 1: + raise RuntimeError( + "Expected to find only one value for skyin " + "and skyout in the APCORR reference file for " + "a given selector." + ) + params["bkg_aperture_inner_radius"] = skyins[0] / pixel_scale_ratio + params["bkg_aperture_outer_radius"] = skyouts[0] / pixel_scale_ratio + + return params + + @lazyproperty + def abvega_offset(self): + """ + Offset to convert from AB to Vega magnitudes. + + The value represents m_AB - m_Vega. + """ + + # TODO (@bmorris3): replace or remove + log.warning( + "ABVEGAOFFSET reference file not yet available. " + "Catalog Vega magnitudes are not correct." + ) + return 0.0 + + # TODO (@bmorris3): replace or remove + # if self.abvegaoffset_filename is None: + # log.warning( + # "ABVEGAOFFSET reference file was not input. " + # "Catalog Vega magnitudes are not correct." + # ) + # return 0.0 + # + # if self.instrument in ("NIRCAM", "NIRISS"): + # selector = {"filter": self.filtername, "pupil": self.pupil} + # elif self.instrument == "MIRI": + # selector = {"filter": self.filtername} + # elif self.instrument == "FGS": + # selector = {"detector": self.detector} + # else: + # raise RuntimeError(f"{self.instrument} is not a valid instrument") + # + # # TODO (@bmorris3): replace or remove + # # abvegaoffset_model = ABVegaOffsetModel(self.abvegaoffset_filename) + # # offsets_table = abvegaoffset_model.abvega_offset + # + # try: + # mask_idx = [offsets_table[key] == value for key, value in selector.items()] + # except KeyError as badkey: + # raise KeyError( + # f"{badkey} not found in ABVEGAOFFSET reference " + # f"file {self.abvegaoffset_filename}" + # ) + # + # row = offsets_table[np.logical_and.reduce(mask_idx)] + # + # if len(row) == 0: + # raise RuntimeError( + # "Did not find matching row in ABVEGAOFFSET " + # f"reference file {self.abvegaoffset_filename}" + # ) + # if len(row) > 1: + # raise RuntimeError( + # "Found more than one matching row in " + # "ABVEGAOFFSET reference file " + # f"{self.abvegaoffset_filename}" + # ) + # + # abvega_offset = row["abvega_offset"][0] + # log.info(f"AB to Vega magnitude offset {abvega_offset:.5f}") + # abvegaoffset_model.close() + # return abvega_offset diff --git a/romancal/source_catalog/source_catalog.py b/romancal/source_catalog/source_catalog.py new file mode 100644 index 000000000..4d9f93bac --- /dev/null +++ b/romancal/source_catalog/source_catalog.py @@ -0,0 +1,1036 @@ +""" +Module to calculate the source catalog. +""" + +import logging +import warnings + +import astropy.units as u +import numpy as np +from astropy.convolution import Gaussian2DKernel +from astropy.nddata.utils import NoOverlapError, extract_array +from astropy.stats import SigmaClip, gaussian_fwhm_to_sigma +from astropy.table import QTable +from astropy.utils import lazyproperty +from astropy.utils.exceptions import AstropyUserWarning +from photutils.aperture import CircularAnnulus, CircularAperture, aperture_photometry +from photutils.segmentation import SourceCatalog +from roman_datamodels.datamodels import ImageModel +from scipy import ndimage +from scipy.spatial import KDTree + +from ._wcs_helpers import pixel_scale_angle_at_skycoord + +log = logging.getLogger(__name__) +log.setLevel(logging.DEBUG) + + +class RomanSourceCatalog: + """ + Class for the Roman source catalog. + + Parameters + ---------- + model : `ImageModel` + The input `ImageModel`. The data is assumed to be + background subtracted. + + segment_image : `~photutils.segmentation.SegmentationImage` + A 2D segmentation image, with the same shape as the input data, + where sources are marked by different positive integer values. + A value of zero is reserved for the background. + + convolved_data : data : 2D `~numpy.ndarray` + The 2D array used to calculate the source centroid and + morphological properties. + + kernel_fwhm : float + The full-width at half-maximum (FWHM) of the 2D Gaussian kernel. + This is needed to calculate the DAOFind sharpness and roundness + properties (DAOFind uses a special kernel that sums to zero). + + aperture_params : `dict` + A dictionary containing the aperture parameters (radii, aperture + corrections, and background annulus inner and outer radii). + + abvega_offset : float + Offset to convert from AB to Vega magnitudes. The value + represents m_AB - m_Vega. + + ci_star_thresholds : array-like of 2 floats + The concentration index thresholds for determining whether + a source is a star. The first threshold corresponds to the + concentration index calculated from the smallest and middle + aperture radii (see ``aperture_params``). The second threshold + corresponds to the concentration index calculated from the + middle and largest aperture radii. An object is considered + extended if both concentration indices are greater than the + corresponding thresholds, otherwise it is considered a star. + + Notes + ----- + ``model.err`` is assumed to be the total error array corresponding + to the input science ``model.data`` array. It is assumed to include + *all* sources of error, including the Poisson error of the sources, + and have the same shape and units as the science data array. + """ + + def __init__( + self, + model, + segment_img, + convolved_data, + kernel_fwhm, + aperture_params, + abvega_offset, + ci_star_thresholds, + ): + + if not isinstance(model, ImageModel): + raise ValueError("The input model must be a ImageModel.") + self.model = model # background was previously subtracted + + self.segment_img = segment_img + self.convolved_data = convolved_data + self.kernel_sigma = kernel_fwhm * gaussian_fwhm_to_sigma + self.aperture_params = aperture_params + self.abvega_offset = abvega_offset + + if len(ci_star_thresholds) != 2: + raise ValueError("ci_star_thresholds must contain only 2 items") + self.ci_star_thresholds = ci_star_thresholds + + self.n_sources = len(self.segment_img.labels) + self.aperture_ee = self.aperture_params["aperture_ee"] + self.n_aper = len(self.aperture_ee) + self.wcs = self.model.meta.wcs + self.column_desc = {} + self._xpeak = None + self._ypeak = None + self.meta = {} + + def convert_to_jy(self): + """ + Convert the data and errors from MJy/sr to Jy and convert to + `~astropy.unit.Quantity` objects. + """ + # TODO (@bmorris3): replace or remove + # in_unit = 'MJy/sr' + # if (self.model.meta.bunit_data != in_unit or + # self.model.meta.bunit_err != in_unit): + # raise ValueError('data and err are expected to be in units of ' + # 'MJy/sr') + # unit = u.Jy + # self.model.meta['photometry'] = dict(pixelarea_steradians=0.11 * u.arcsec**2) + # to_jy = 1.e6 * self.model.meta.photometry.pixelarea_steradians + # self.model.data = self.model.data + # self.model.err *= to_jy + # self.model.data <<= unit + # self.model.err <<= unit + # self.model.meta.bunit_data = unit.name + # self.model.meta.bunit_err = unit.name + pass + + def convert_from_jy(self): + """ + Convert the data and errors from Jy to MJy/sr and change from + `~astropy.unit.Quantity` objects to `~numpy.ndarray`. + """ + # TODO (@bmorris3): replace or remove + # to_mjy_sr = 1.e6 * self.model.meta.photometry.pixelarea_steradians + # self.model.data /= to_mjy_sr + # self.model.err /= to_mjy_sr + # + # self.model.data = self.model.data.value # remove units + # self.model.err = self.model.err.value # remove units + # + # self.model.meta.bunit_data = 'MJy/sr' + # self.model.meta.bunit_err = 'MJy/sr' + + pass + + @staticmethod + def convert_flux_to_abmag(flux, flux_err): + """ + Convert flux (and error) to AB magnitude (and error). + + Parameters + ---------- + flux, flux_err : `~astropy.unit.Quantity` + The input flux and error arrays in units of Jy. + + Returns + ------- + abmag, abmag_err : `~astropy.ndarray` + The output AB magnitude and error arrays. + """ + # ignore RunTimeWarning if flux or flux_err contains NaNs + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=RuntimeWarning) + + abmag = -2.5 * np.log10(flux.value) + 8.9 + abmag_err = 2.5 * np.log10(1.0 + (flux_err.value / flux.value)) + + # handle negative fluxes + idx = flux.value < 0 + abmag[idx] = np.nan + abmag_err[idx] = np.nan + + return abmag, abmag_err + + @lazyproperty + def segment_colnames(self): + """ + A dictionary of the output table column names and descriptions + for the segment catalog. + """ + desc = {} + desc["label"] = "Unique source identification label number" + desc["xcentroid"] = "X pixel value of the source centroid (0 indexed)" + desc["ycentroid"] = "Y pixel value of the source centroid (0 indexed)" + desc["sky_centroid"] = "ICRS Sky coordinate of the source centroid" + desc["isophotal_flux"] = "Isophotal flux" + desc["isophotal_flux_err"] = "Isophotal flux error" + desc["isophotal_abmag"] = "Isophotal AB magnitude" + desc["isophotal_abmag_err"] = "Isophotal AB magnitude error" + desc["isophotal_vegamag"] = "Isophotal Vega magnitude" + desc["isophotal_vegamag_err"] = "Isophotal Vega magnitude error" + desc["isophotal_area"] = "Isophotal area" + desc["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" + ) + desc["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" + ) + desc["ellipticity"] = ( + "1 minus the ratio of the 1-sigma lengths of " + "the semimajor and semiminor axes" + ) + desc["orientation"] = ( + "The angle (degrees) between the positive X " + "axis and the major axis (increases " + "counter-clockwise)" + ) + desc["sky_orientation"] = ( + "The position angle (degrees) from North " "of the major axis" + ) + desc["sky_bbox_ll"] = ( + "Sky coordinate of the lower-left vertex of " + "the minimal bounding box of the source" + ) + desc["sky_bbox_ul"] = ( + "Sky coordinate of the upper-left vertex of " + "the minimal bounding box of the source" + ) + desc["sky_bbox_lr"] = ( + "Sky coordinate of the lower-right vertex of " + "the minimal bounding box of the source" + ) + desc["sky_bbox_ur"] = ( + "Sky coordinate of the upper-right vertex of " + "the minimal bounding box of the source" + ) + + self.column_desc.update(desc) + + return list(desc.keys()) + + def set_segment_properties(self): + """ + Calculate the segment-based source photometry and morphologies. + + The values are set as dynamic attributes. + """ + segm_cat = SourceCatalog( + self.model.data, + self.segment_img, + convolved_data=self.convolved_data << u.electron / u.s, + error=self.model.err, + wcs=self.wcs, + ) + self._xpeak = segm_cat.maxval_xindex + self._ypeak = segm_cat.maxval_yindex + + self.meta.update(segm_cat.meta) + for key in ("sklearn", "matplotlib"): + self.meta["version"].pop(key) + + # rename some columns in the output catalog + prop_names = {} + prop_names["isophotal_flux"] = "segment_flux" + prop_names["isophotal_flux_err"] = "segment_fluxerr" + prop_names["isophotal_area"] = "area" + + for column in self.segment_colnames: + # define the property name + prop_name = prop_names.get(column, column) + try: + value = getattr(segm_cat, prop_name) + except AttributeError: + value = getattr(self, prop_name) + setattr(self, column, value) + + @lazyproperty + def xypos(self): + """ + The (x, y) source positions, defined from the segmentation + image. + """ + return np.transpose((self.xcentroid, self.ycentroid)) + + @lazyproperty + def _xypos_finite(self): + """ + The (x, y) source positions, where non-finite positions are + set to a large negative value. + + At this position the aperture will not overlap the data, thus + returning NaN fluxes and errors. + """ + xypos = self.xypos.copy() + nanmask = ~np.isfinite(xypos) + xypos[nanmask] = -1000.0 + return xypos + + @lazyproperty + def _xypos_nonfinite_mask(self): + """ + A 1D boolean mask where `True` values denote sources where + either the xcentroid or the ycentroid is not finite. + """ + return ~np.isfinite(self.xypos).all(axis=1) + + @lazyproperty + def _isophotal_abmag(self): + """ + The isophotal AB magnitude and error. + """ + return self.convert_flux_to_abmag(self.isophotal_flux, self.isophotal_flux_err) + + @lazyproperty + def isophotal_abmag(self): + """ + The isophotal AB magnitude. + """ + return self._isophotal_abmag[0] + + @lazyproperty + def isophotal_abmag_err(self): + """ + The isophotal AB magnitude error. + """ + return self._isophotal_abmag[1] + + @lazyproperty + def isophotal_vegamag(self): + """ + The isophotal Vega magnitude. + """ + return self.isophotal_abmag - self.abvega_offset + + @lazyproperty + def isophotal_vegamag_err(self): + """ + The isophotal Vega magnitude error. + """ + return self.isophotal_abmag_err + + @lazyproperty + def sky_orientation(self): + """ + The orientation of the source major axis as the position angle + in degrees measured East of North. + """ + skycoord = self.wcs.pixel_to_world(0 * u.pix, 0 * u.pix) + _, _, angle = pixel_scale_angle_at_skycoord(skycoord, self.wcs) + + return (180.0 * u.deg) - angle + self.orientation + + def _make_aperture_colnames(self, name): + """ + Make the aperture column names. + + There are separate columns for the flux/magnitude for each of + the encircled energies and the total flux/magnitude. + + Parameters + ---------- + name : {'flux', 'abmag', 'vegamag'} + The name type of the column. + + Returns + ------- + colnames : list of str + A list of the output column names. + """ + colnames = [] + for aper_ee in self.aperture_ee: + basename = f"aper{aper_ee}_{name}" + colnames.append(basename) + colnames.append(f"{basename}_err") + colnames.extend([f"aper_total_{name}", f"aper_total_{name}_err"]) + + return colnames + + def _make_aperture_descriptions(self, name): + """ + Make aperture column descriptions. + + Parameters + ---------- + name : {'flux', 'abmag', 'vegamag'} + The name type of the column. + + Returns + ------- + descriptions : list of str + A list of the output column descriptions. + """ + if name == "flux": + ftype = "Flux" + ftype2 = "flux" + elif name == "abmag": + ftype = ftype2 = "AB magnitude" + elif name == "vegamag": + ftype = ftype2 = "Vega magnitude" + + desc = [] + for aper_ee in self.aperture_ee: + desc.append( + f"{ftype} within the {aper_ee}% encircled energy " "circular aperture" + ) + desc.append( + f"{ftype} error within the {aper_ee}% encircled " + "energy circular aperture" + ) + + desc.append( + f"Total aperture-corrected {ftype2} based on the " + f"{self.aperture_ee[-1]}% encircled energy circular " + "aperture; should be used only for unresolved sources." + ) + desc.append( + f"Total aperture-corrected {ftype2} error based on the " + f"{self.aperture_ee[-1]}% encircled energy circular " + "aperture; should be used only for unresolved sources." + ) + + return desc + + @lazyproperty + def aperture_flux_colnames(self): + """ + The aperture flux column names. + """ + return self._make_aperture_colnames("flux") + + @lazyproperty + def aperture_flux_descriptions(self): + """ + The aperture flux column descriptions. + """ + return self._make_aperture_descriptions("flux") + + @lazyproperty + def aperture_abmag_colnames(self): + """ + The aperture AB magnitude column names. + """ + return self._make_aperture_colnames("abmag") + + @lazyproperty + def aperture_abmag_descriptions(self): + """ + The aperture AB magnitude column descriptions. + """ + return self._make_aperture_descriptions("abmag") + + @lazyproperty + def aperture_vegamag_colnames(self): + """ + The aperture Vega magnitude column names. + """ + return self._make_aperture_colnames("vegamag") + + @lazyproperty + def aperture_vegamag_descriptions(self): + """ + The aperture Vega magnitude column descriptions. + """ + return self._make_aperture_descriptions("vegamag") + + @lazyproperty + def aperture_colnames(self): + """ + A dictionary of the output table column names and descriptions + for the aperture catalog. + """ + desc = {} + desc["aper_bkg_flux"] = ( + "The local background value calculated as " + "the sigma-clipped median value in the " + "background annulus aperture" + ) + desc["aper_bkg_flux_err"] = ( + "The standard error of the " "sigma-clipped median background value" + ) + + for idx, colname in enumerate(self.aperture_flux_colnames): + desc[colname] = self.aperture_flux_descriptions[idx] + for idx, colname in enumerate(self.aperture_abmag_colnames): + desc[colname] = self.aperture_abmag_descriptions[idx] + for idx, colname in enumerate(self.aperture_vegamag_colnames): + desc[colname] = self.aperture_vegamag_descriptions[idx] + + self.column_desc.update(desc) + + return list(desc.keys()) + + @lazyproperty + def _aper_local_background(self): + """ + Estimate the local background and error using a circular annulus + aperture. + + The local background is the sigma-clipped median value in the + annulus. The background error is the standard error of the + median, sqrt(pi / 2N) * std. + """ + bkg_aper = CircularAnnulus( + self._xypos_finite, + self.aperture_params["bkg_aperture_inner_radius"], + self.aperture_params["bkg_aperture_outer_radius"], + ) + bkg_aper_masks = bkg_aper.to_mask(method="center") + sigclip = SigmaClip(sigma=3.0) + + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=RuntimeWarning) + warnings.simplefilter("ignore", category=AstropyUserWarning) + + nvalues = [] + bkg_median = [] + bkg_std = [] + for mask in bkg_aper_masks: + bkg_data = mask.get_values(self.model.data.value) + values = sigclip(bkg_data, masked=False) + nvalues.append(values.size) + bkg_median.append(np.median(values)) + bkg_std.append(np.std(values)) + + nvalues = np.array(nvalues) + bkg_median = np.array(bkg_median) + # standard error of the median + bkg_median_err = np.sqrt(np.pi / (2.0 * nvalues)) * np.array(bkg_std) + + bkg_median <<= self.model.data.unit + bkg_median_err <<= self.model.data.unit + + return bkg_median, bkg_median_err + + @lazyproperty + def aper_bkg_flux(self): + """ + The aperture local background flux (per pixel). + """ + return self._aper_local_background[0] + + @lazyproperty + def aper_bkg_flux_err(self): + """ + The aperture local background flux error (per pixel). + """ + return self._aper_local_background[1] + + def set_aperture_properties(self): + """ + Calculate the aperture photometry. + + The values are set as dynamic attributes. + """ + apertures = [ + CircularAperture(self._xypos_finite, radius) + for radius in self.aperture_params["aperture_radii"] + ] + aper_phot = aperture_photometry( + self.model.data, apertures, error=self.model.err + ) + + for i, aperture in enumerate(apertures): + flux_col = f"aperture_sum_{i}" + flux_err_col = f"aperture_sum_err_{i}" + + # subtract the local background measured in the annulus + aper_phot[flux_col] -= self.aper_bkg_flux * aperture.area + + flux = aper_phot[flux_col] + flux_err = aper_phot[flux_err_col] + abmag, abmag_err = self.convert_flux_to_abmag(flux, flux_err) + vegamag = abmag - self.abvega_offset + vegamag_err = abmag_err + + idx0 = 2 * i + idx1 = (2 * i) + 1 + setattr(self, self.aperture_flux_colnames[idx0], flux) + setattr(self, self.aperture_flux_colnames[idx1], flux_err) + setattr(self, self.aperture_abmag_colnames[idx0], abmag) + setattr(self, self.aperture_abmag_colnames[idx1], abmag_err) + setattr(self, self.aperture_vegamag_colnames[idx0], vegamag) + setattr(self, self.aperture_vegamag_colnames[idx1], vegamag_err) + + @lazyproperty + def extras_colnames(self): + """ + A dictionary of the output table column names and descriptions + for the additional catalog values. + """ + desc = {} + for idx, colname in enumerate(self.ci_colnames): + desc[colname] = self.ci_colname_descriptions[idx] + + desc["is_extended"] = "Flag indicating whether the source is extended" + desc["sharpness"] = "The DAOFind source sharpness statistic" + desc["roundness"] = "The DAOFind source roundness statistic" + desc["nn_label"] = "The label number of the nearest neighbor" + desc["nn_dist"] = "The distance in pixels to the nearest neighbor" + self.column_desc.update(desc) + + return list(desc.keys()) + + @lazyproperty + def _ci_ee_indices(self): + """ + The EE indices for the concentration indices. + """ + # NOTE: the EE values are always in increasing order + return (0, 1), (1, 2), (0, 2) + + @lazyproperty + def ci_colnames(self): + """ + The column names of the three concentration indices. + """ + return [ + f"CI_{self.aperture_ee[j]}_{self.aperture_ee[i]}" + for (i, j) in self._ci_ee_indices + ] + + @lazyproperty + def ci_colname_descriptions(self): + """ + The concentration indices column descriptions. + """ + return [ + "Concentration index calculated as " + f"({self.aperture_flux_colnames[2*j]} / " + f"{self.aperture_flux_colnames[2*i]})" + for (i, j) in self._ci_ee_indices + ] + + @lazyproperty + def concentration_indices(self): + """ + A list of concentration indices, calculated as the flux + ratios of: + + * the middle / smallest aperture radii/EE, + e.g., CI_50_30 = aper50_flux / aper30_flux + * the largest / middle aperture radii/EE, + e.g., CI_70_50 = aper70_flux / aper50_flux + * the largest / smallest aperture radii/EE, + e.g., CI_70_30 = aper70_flux / aper30_flux + """ + fluxes = [ + (self.aperture_flux_colnames[2 * j], self.aperture_flux_colnames[2 * i]) + for (i, j) in self._ci_ee_indices + ] + return [ + getattr(self, flux1).value / getattr(self, flux2).value + for flux1, flux2 in fluxes + ] + + def set_ci_properties(self): + """ + Set the concentration indices as dynamic attributes. + """ + for name, value in zip(self.ci_colnames, self.concentration_indices): + setattr(self, name, value) + + @lazyproperty + def is_extended(self): + """ + Boolean indicating whether the source is extended. + """ + mask1 = self.concentration_indices[0] > self.ci_star_thresholds[0] + mask2 = self.concentration_indices[1] > self.ci_star_thresholds[1] + return np.logical_and(mask1, mask2) + + @lazyproperty + def _kernel_size(self): + """ + The DAOFind kernel size (in both x and y dimensions). + """ + # always odd + return 2 * int(max(2.0, 1.5 * self.kernel_sigma)) + 1 + + @lazyproperty + def _kernel_center(self): + """ + The DAOFind kernel x/y center. + """ + return (self._kernel_size - 1) // 2 + + @lazyproperty + def _kernel_mask(self): + """ + The DAOFind kernel circular mask. + NOTE: 1=good pixels, 0=masked pixels + """ + yy, xx = np.mgrid[0 : self._kernel_size, 0 : self._kernel_size] + radius = np.sqrt( + (xx - self._kernel_center) ** 2 + (yy - self._kernel_center) ** 2 + ) + return (radius <= max(2.0, 1.5 * self.kernel_sigma)).astype(int) + + @lazyproperty + def _daofind_kernel(self): + """ + The DAOFind kernel, a 2D circular Gaussian normalized to have + zero sum. + """ + size = self._kernel_size + kernel = Gaussian2DKernel(self.kernel_sigma, x_size=size, y_size=size).array + kernel /= np.max(kernel) + kernel *= self._kernel_mask + + # normalize the kernel to zero sum + npixels = self._kernel_mask.sum() + denom = np.sum(kernel**2) - (np.sum(kernel) ** 2 / npixels) + return ((kernel - (kernel.sum() / npixels)) / denom) * self._kernel_mask + + @lazyproperty + def _daofind_convolved_data(self): + """ + The DAOFind convolved data. + """ + return ndimage.convolve( + self.model.data.value, self._daofind_kernel, mode="constant", cval=0.0 + ) + + @lazyproperty + def _daofind_cutout(self): + """ + 3D array containing 2D cutouts centered on each source from the + input data. + + The cutout size always matches the size of the DAOFind kernel, + which has odd dimensions. + """ + cutout = [] + for xcen, ycen in zip(*np.transpose(self._xypos_finite)): + try: + cutout_ = extract_array( + self.model.data, + self._daofind_kernel.shape, + (ycen, xcen), + fill_value=0.0, + ) + except NoOverlapError: + cutout_ = np.zeros(self._daofind_kernel.shape) + cutout.append(cutout_) + + return np.array(cutout) # all cutouts are the same size + + @lazyproperty + def _daofind_cutout_conv(self): + """ + 3D array containing 2D cutouts centered on each source from the + DAOFind convolved data. + + The cutout size always matches the size of the DAOFind kernel, + which has odd dimensions. + """ + cutout = [] + for xcen, ycen in zip(*np.transpose(self._xypos_finite)): + try: + cutout_ = extract_array( + self._daofind_convolved_data, + self._daofind_kernel.shape, + (ycen, xcen), + fill_value=0.0, + ) + except NoOverlapError: + cutout_ = np.zeros(self._daofind_kernel.shape) + cutout.append(cutout_) + + return np.array(cutout) # all cutouts are the same size + + @lazyproperty + def sharpness(self): + """ + The DAOFind source sharpness statistic. + + The sharpness statistic measures the ratio of the difference + between the height of the central pixel and the mean of the + surrounding non-bad pixels to the height of the best fitting + Gaussian function at that point. + + Stars generally have a ``sharpness`` between 0.2 and 1.0. + """ + npixels = self._kernel_mask.sum() - 1 # exclude the peak pixel + data_masked = self._daofind_cutout * self._kernel_mask + data_peak = self._daofind_cutout[:, self._kernel_center, self._kernel_center] + conv_peak = self._daofind_cutout_conv[ + :, self._kernel_center, self._kernel_center + ] + + data_mean = (np.sum(data_masked, axis=(1, 2)) - data_peak) / npixels + + with warnings.catch_warnings(): + # ignore 0 / 0 for non-finite xypos + warnings.simplefilter("ignore", category=RuntimeWarning) + return (data_peak - data_mean) / conv_peak + + @lazyproperty + def roundness(self): + """ + The DAOFind source roundness statistic based on symmetry. + + The roundness characteristic computes the ratio of a measure of + the bilateral symmetry of the object to a measure of the + four-fold symmetry of the object. + + "Round" objects have a ``roundness`` close to 0, generally + between -1 and 1. + """ + # set the central (peak) pixel to zero + cutout = self._daofind_cutout_conv.copy() + cutout[:, self._kernel_center, self._kernel_center] = 0.0 + + # calculate the four roundness quadrants + quad1 = cutout[:, 0 : self._kernel_center + 1, self._kernel_center + 1 :] + quad2 = cutout[:, 0 : self._kernel_center, 0 : self._kernel_center + 1] + quad3 = cutout[:, self._kernel_center :, 0 : self._kernel_center] + quad4 = cutout[:, self._kernel_center + 1 :, self._kernel_center :] + + axis = (1, 2) + sum2 = ( + -quad1.sum(axis=axis) + + quad2.sum(axis=axis) + - quad3.sum(axis=axis) + + quad4.sum(axis=axis) + ) + sum2[sum2 == 0] = 0.0 + + sum4 = np.abs(cutout).sum(axis=axis) + sum4[sum4 == 0] = np.nan + + with warnings.catch_warnings(): + # ignore 0 / 0 for non-finite xypos + warnings.simplefilter("ignore", category=RuntimeWarning) + return 2.0 * sum2 / sum4 + + @lazyproperty + def _kdtree_query(self): + """ + The distance in pixels to the nearest neighbor and its index. + """ + if self.n_sources == 1: + return [np.nan], [np.nan] + + # non-finite xypos causes memory errors on linux, but not MacOS + tree = KDTree(self._xypos_finite) + qdist, qidx = tree.query(self._xypos_finite, k=[2]) + return np.transpose(qdist)[0], np.transpose(qidx)[0] + + @lazyproperty + def nn_label(self): + """ + The label number of the nearest neighbor. + + A label value of -1 is returned if there is only one detected + source and for sources with a non-finite xcentroid or ycentroid. + """ + if self.n_sources == 1: + return -1 + + nn_label = self.label[self._kdtree_query[1]] + # assign a label of -1 for non-finite xypos + nn_label[self._xypos_nonfinite_mask] = -1 + + return nn_label + + @lazyproperty + def nn_dist(self): + """ + The distance in pixels to the nearest neighbor. + """ + nn_dist = self._kdtree_query[0] + if self.n_sources == 1: + # NaN if only one detected source + return nn_dist * u.pixel + + # assign a distance of np.nan for non-finite xypos + nn_dist[self._xypos_nonfinite_mask] = np.nan + return nn_dist * u.pixel + + @lazyproperty + def aper_total_flux(self): + """ + The aperture-corrected total flux for sources, based on the flux + in largest aperture. + + The aperture-corrected total flux should be used only for + unresolved sources. + """ + idx = self.n_aper - 1 # apcorr for the largest EE (largest radius) + flux = self.aperture_params["aperture_corrections"][idx] * getattr( + self, self.aperture_flux_colnames[idx * 2] + ) + return flux + + @lazyproperty + def aper_total_flux_err(self): + """ + The aperture-corrected total flux error for sources, + based on the flux in largest aperture. + + The aperture-corrected total flux error should be used only for + unresolved sources. + """ + idx = self.n_aper - 1 # apcorr for the largest EE (largest radius) + flux_err = self.aperture_params["aperture_corrections"][idx] * getattr( + self, self.aperture_flux_colnames[idx * 2 + 1] + ) + return flux_err + + @lazyproperty + def _abmag_total(self): + """ + The total AB magnitude and error. + """ + return self.convert_flux_to_abmag( + self.aper_total_flux, self.aper_total_flux_err + ) + + @lazyproperty + def aper_total_abmag(self): + """ + The aperture-corrected total AB magnitude. + + The aperture-corrected total magnitude should be used only for + unresolved sources. + """ + return self._abmag_total[0] + + @lazyproperty + def aper_total_abmag_err(self): + """ + The aperture-corrected total AB magnitude error. + + The aperture-corrected total magnitude error should be used only + for unresolved sources. + """ + return self._abmag_total[1] + + @lazyproperty + def aper_total_vegamag(self): + """ + The aperture-corrected total Vega magnitude. + + The aperture-corrected total magnitude should be used only for + unresolved sources. + """ + return self.aper_total_abmag - self.abvega_offset + + @lazyproperty + def aper_total_vegamag_err(self): + """ + The aperture-corrected total Vega magnitude error. + + The aperture-corrected total magnitude error should be used only + for unresolved sources. + """ + return self.aper_total_abmag_err + + @lazyproperty + def colnames(self): + """ + The column name order for the final source catalog. + """ + colnames = self.segment_colnames[0:4] + colnames.extend(self.aperture_colnames) + colnames.extend(self.extras_colnames) + colnames.extend(self.segment_colnames[4:]) + return colnames + + @staticmethod + def format_columns(catalog): + """ + Format the values in the output catalog. + + Parameters + ---------- + catalog : `~astropy.table.Table` + The catalog to format. + + Returns + ------- + result : `~astropy.table.Table` + The formated catalog. + """ + # output formatting requested by the JPWG (2020.02.05) + for colname in catalog.colnames: + if colname in ("xcentroid", "ycentroid") or "CI_" in colname: + catalog[colname].info.format = ".4f" + if "flux" in colname: + catalog[colname].info.format = ".6e" + if ( + "abmag" in colname + or "vegamag" in colname + or colname in ("nn_dist", "sharpness", "roundness") + ): + catalog[colname].info.format = ".6f" + if colname in ( + "semimajor_sigma", + "semiminor_sigma", + "ellipticity", + "orientation", + "sky_orientation", + ): + catalog[colname].info.format = ".6f" + + return catalog + + @lazyproperty + def catalog(self): + """ + The final source catalog. + """ + self.convert_to_jy() + self.set_segment_properties() + self.set_aperture_properties() + self.set_ci_properties() + + catalog = QTable() + for column in self.colnames: + catalog[column] = getattr(self, column) + catalog[column].info.description = self.column_desc[column] + + catalog = self.format_columns(catalog) + + # update metadata + self.meta["aperture_params"] = self.aperture_params + self.meta["abvega_offset"] = self.abvega_offset + catalog.meta.update(self.meta) + + # reset units on input model back to MJy/sr + self.convert_from_jy() + + return catalog diff --git a/romancal/source_catalog/source_catalog_step.py b/romancal/source_catalog/source_catalog_step.py new file mode 100755 index 000000000..40e9af1fd --- /dev/null +++ b/romancal/source_catalog/source_catalog_step.py @@ -0,0 +1,123 @@ +""" +Module for the source catalog step. +""" + +import os + +import numpy as np +from crds.core.exceptions import CrdsLookupError +from roman_datamodels import datamodels + +from romancal.stpipe import RomanStep + +from .detection import RomanBackground, RomanSourceFinder, convolve_data +from .reference_data import ReferenceData +from .source_catalog import RomanSourceCatalog + +__all__ = ["SourceCatalogStep"] + + +class SourceCatalogStep(RomanStep): + """ + Create a final catalog of source photometry and morphologies. + + Parameters + ----------- + input : str or `ImageModel` + Path to an ASDF file, or an `ImageModel`. + """ + + class_alias = "source_catalog" + + spec = """ + bkg_boxsize = integer(default=1000) # background mesh box size in pixels + kernel_fwhm = float(default=2.0) # Gaussian kernel FWHM in pixels + snr_threshold = float(default=3.0) # SNR threshold above the bkg + npixels = integer(default=25) # min number of pixels in source + deblend = boolean(default=False) # deblend sources? + aperture_ee1 = integer(default=30) # aperture encircled energy 1 + aperture_ee2 = integer(default=50) # aperture encircled energy 2 + aperture_ee3 = integer(default=70) # aperture encircled energy 3 + ci1_star_threshold = float(default=2.0) # CI 1 star threshold + ci2_star_threshold = float(default=1.8) # CI 2 star threshold + suffix = string(default='cat') # Default suffix for output files + """ + + def _get_reffile_paths(self, model): + filepaths = [] + for reffile_type in self.reference_file_types: + try: + filepath = self.get_reference_file(model, reffile_type) + self.log.info( + f"Using {reffile_type.upper()} reference file: " f"{filepath}" + ) + except CrdsLookupError as err: + msg = f"{err} Source catalog will not be created." + self.log.warning(msg) + return None + + filepaths.append(filepath) + return filepaths + + def process(self, input_model): + with datamodels.open(input_model) as model: + reffile_paths = self._get_reffile_paths(model) + aperture_ee = (self.aperture_ee1, self.aperture_ee2, self.aperture_ee3) + + try: + refdata = ReferenceData(model, reffile_paths, aperture_ee) + aperture_params = refdata.aperture_params + + # TODO (@bmorris3): replace or remove + abvega_offset = refdata.abvega_offset + except RuntimeError as err: + msg = f"{err} Source catalog will not be created." + self.log.warning(msg) + return None + + coverage_mask = np.isnan(model.err) # | (model.wht == 0) + bkg = RomanBackground( + model.data, box_size=self.bkg_boxsize, coverage_mask=coverage_mask + ) + model.data -= bkg.background + + threshold = self.snr_threshold * bkg.background_rms + finder = RomanSourceFinder(threshold, self.npixels, deblend=self.deblend) + + convolved_data = convolve_data( + model.data, self.kernel_fwhm, mask=coverage_mask + ) + segment_img = finder(convolved_data, mask=coverage_mask) + if segment_img is None: + return None + + ci_star_thresholds = (self.ci1_star_threshold, self.ci2_star_threshold) + catobj = RomanSourceCatalog( + model, + segment_img, + convolved_data, + self.kernel_fwhm, + aperture_params, + abvega_offset, + ci_star_thresholds, + ) + catalog = catobj.catalog + + # add back background to data so input model is unchanged + model.data += bkg.background + + if self.save_results: + cat_filepath = self.make_output_path(ext=".ecsv") + catalog.write(cat_filepath, format="ascii.ecsv", overwrite=True) + model.meta.source_catalog = os.path.basename(cat_filepath) + self.log.info(f"Wrote source catalog: {cat_filepath}") + + segm_model = datamodels.SegmentationMapModel(segment_img.data) + segm_model.update(model, only="PRIMARY") + segm_model.meta.wcs = model.meta.wcs + segm_model.meta.wcsinfo = model.meta.wcsinfo + self.save_model(segm_model, suffix="segm") + model.meta.segmentation_map = segm_model.meta.filename + self.log.info("Wrote segmentation map: " f"{segm_model.meta.filename}") + + return catalog diff --git a/romancal/source_catalog/tests/__init__.py b/romancal/source_catalog/tests/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/romancal/source_catalog/tests/test_source_catalog.py b/romancal/source_catalog/tests/test_source_catalog.py new file mode 100644 index 000000000..f785219fe --- /dev/null +++ b/romancal/source_catalog/tests/test_source_catalog.py @@ -0,0 +1,159 @@ +import astropy.units as u +import numpy as np +import pytest +from numpy.testing import assert_allclose +from roman_datamodels.datamodels import ImageModel +from roman_datamodels.maker_utils import mk_level2_image + +from romancal.lib import dqflags +from romancal.source_catalog.source_catalog_step import SourceCatalogStep + +DO_NOT_USE = dqflags.pixel["DO_NOT_USE"] +SATURATED = dqflags.pixel["SATURATED"] + + +def mk_image_model( + rate_mean=0, + rate_std=1e-4, + image_shape=(100, 100), + rng=np.random.default_rng(619), +): + l2 = mk_level2_image(shape=image_shape) + l2_im = ImageModel(l2) + l2_im.data = u.Quantity( + rng.normal(loc=rate_mean, scale=rate_std, size=l2_im.data.shape).astype( + np.float32 + ), + l2_im.data.unit, + ) + + # fake a background until `rad` implements the schema: + l2_im.meta["background"] = dict(level=None, subtracted=False, method=None) + + l2_im.meta.cal_step["skymatch"] = "INCOMPLETE" + return l2_im + + +@pytest.fixture +def wfi_model(): + units = u.electron / u.s + rng = np.random.default_rng(seed=123) + data = rng.normal(0, 0.5, size=(101, 101)) + data[20:80, 10:20] = 1.4 + data[20:30, 20:45] = 1.4 + data[20:80, 55:65] = 7.2 + data[70:80, 65:87] = 7.2 + data[45:55, 65:87] = 7.2 + data[20:30, 65:87] = 7.2 + data[55:75, 82:92] = 7.2 + data[25:45, 82:92] = 7.2 + + err = np.abs(data) / 10.0 + + model = mk_image_model(0, 0.5, image_shape=(101, 101), rng=rng) + model.data = data << units + model.err = err << units + + # model.meta.bunit_data = 'e/s' + # model.meta.bunit_err = 'e/s' + + return model + + +# TODO (@bmorris3): replace or remove +# @pytest.fixture +# def wfi_model_without_apcorr(): +# rng = np.random.default_rng(seed=123) +# data = rng.normal(0, 0.5, size=(101, 101)) +# data[20:80, 10:20] = 1.4 +# data[20:30, 20:45] = 1.4 +# data[20:80, 55:65] = 7.2 +# data[70:80, 65:87] = 7.2 +# data[45:55, 65:87] = 7.2 +# data[20:30, 65:87] = 7.2 +# data[55:75, 82:92] = 7.2 +# data[25:45, 82:92] = 7.2 +# +# wht = np.ones(data.shape) +# wht[0:10, :] = 0. +# err = np.abs(data) / 10. +# model = ImageModel(data, wht=wht, err=err) +# model.meta.bunit_data = 'MJy/sr' +# model.meta.bunit_err = 'MJy/sr' +# model.meta.photometry.pixelarea_steradians = 1.0 +# model.meta.wcs = make_gwcs(data.shape) +# model.meta.wcsinfo = { +# 'ctype1': 'RA---TAN', +# 'ctype2': 'DEC--TAN', +# 'dec_ref': 11.99875540218638, +# 'ra_ref': 22.02351763251896, +# 'roll_ref': 0.005076934167039675, +# 'v2_ref': 86.039011, +# 'v3_ref': -493.385704, +# 'v3yangle': -0.07385127, +# 'vparity': -1, +# 'wcsaxes': 2, +# 'crpix1': 50, +# 'crpix2': 50} +# model.meta.instrument = { +# 'channel': 'LONG', +# 'detector': 'NRCALONG', +# 'filter': 'F2550WR', +# 'lamp_mode': 'NONE', +# 'module': 'A', +# 'name': 'NIRCAM', +# 'pupil': 'CLEAR'} +# model.meta.exposure.type = 'NRC_IMAGE' +# model.meta.observation.date = '2021-01-01' +# model.meta.observation.time = '00:00:00' +# +# return model + + +@pytest.mark.parametrize("npixels, nsources", ((5, 2), (1000, 1), (5000, 0))) +def test_source_catalog(wfi_model, npixels, nsources): + + step = SourceCatalogStep( + snr_threshold=0.5, + npixels=npixels, + bkg_boxsize=50, + kernel_fwhm=2.0, + save_results=False, + ) + cat = step.run(wfi_model) + if cat is None: + assert nsources == 0 + else: + assert len(cat) == nsources + min_snr = np.min(cat["isophotal_flux"] / cat["isophotal_flux_err"]) + assert min_snr >= 100 + + +# TODO (@bmorris3): replace or remove +# def test_model_without_apcorr_data(wfi_model_without_apcorr): +# step = SourceCatalogStep(save_results=False) +# cat = step.run(wfi_model_without_apcorr) +# assert cat is None + + +def test_input_model_reset(wfi_model): + """Changes to input model data are made in SourceCatalogStep - make sure that + these changes are correctly reversed so the input model data/err arrays + remain unchanged after processing (to avoid copying datamodel), + and that the units are in MJy/sr before and after.""" + + original_data = wfi_model.data.copy() + original_err = wfi_model.err.copy() + + step = SourceCatalogStep( + snr_threshold=0.5, + npixels=5, + bkg_boxsize=50, + kernel_fwhm=2.0, + save_results=False, + ) + + step.run(wfi_model) + + assert_allclose(original_data, wfi_model.data, atol=5.0e-7) + assert_allclose(original_err, wfi_model.err, 5.0e-7)