RCAL-1380: Populate SFD-estimated extinction in catalogs.

pyproject.toml

@@ -25,7 +25,8 @@ dependencies = [
     # the roman_datamodels pin to the released version. If there
     # is a need to change this to main then we need a new
     # roman_datamodels release.
-    "roman_datamodels>=0.31.0, <0.32",
+    # "roman_datamodels>=0.31.0, <0.32",
+    "roman_datamodels@git+https://github.com/spacetelescope/roman_datamodels@main",
     # "romanisim>=0.13.1,<0.14",
     "romanisim@git+https://github.com/spacetelescope/romanisim@main",
     "asdf>=4.1.0,<6",

romancal/multiband_catalog/tests/test_multiband_catalog.py

@@ -229,6 +229,24 @@ def test_multiband_catalog(
     shared_tests(result, cat, library_model, save_results, function_jail)
 
 
+def test_multiband_catalog_populates_dust_ebv(library_model):
+    """Ensure the joined multiband catalog contains the detection-level dust_ebv."""
+    step = MultibandCatalogStep()
+    result = step.call(
+        library_model,
+        bkg_boxsize=50,
+        snr_threshold=3,
+        npixels=10,
+        fit_psf=False,
+        save_results=False,
+        deblend=True,
+    )
+    cat = result.source_catalog
+    assert "dust_ebv" in cat.colnames
+    assert len(cat["dust_ebv"]) == len(cat)
+    assert cat["dust_ebv"].dtype == np.float32
+
+
 @pytest.mark.parametrize("save_results", (True, False))
 def test_multiband_catalog_no_detections(library_model, save_results, function_jail):
     step = MultibandCatalogStep()

romancal/regtest/test_multiband_catalog.py

@@ -35,6 +35,7 @@
     "aper02_f213m_flux",  # DMS539 PSF-matched photometry
     "segment_f213m_flux",  # DMS539
     "kron_f213m_flux",  # DMS539
+    "dust_ebv",  # dust extinction values
 ]
 
 

romancal/regtest/test_source_catalog.py

@@ -63,6 +63,7 @@ def fields(catalog):
         "segment_flux_err",  # DMS386 flux uncertainties
         "is_extended",  # DMS376 type of source
         "warning_flags",  # DMS387 dq_flags
+        "dust_ebv",
     ),
 )
 def test_has_field(fields, field):

romancal/source_catalog/_column_schema.py

@@ -206,6 +206,9 @@ def column_names(self):
             "psf_flags",
             "psf_gof",
         ]
+        dust_colnames = [
+            "dust_ebv",
+        ]
 
         det_colnames = []
         det_colnames.extend(segm_colnames)
@@ -230,6 +233,7 @@ def column_names(self):
             colnames.extend(flag_columns)
             if self.fit_psf:
                 colnames.extend(psf_flags_colnames)
+            colnames.extend(dust_colnames)
 
         elif self.cat_type == "forced_det":
             colnames = ["label"]  # needed to join the forced catalogs
@@ -248,6 +252,7 @@ def column_names(self):
             colnames.extend(skywin_colnames)
             colnames.extend(det_colnames)
             colnames.extend(flag_columns)
+            colnames.extend(dust_colnames)
 
         elif self.cat_type == "dr_band":
             colnames = self.band_colnames.copy()

romancal/source_catalog/_source_catalog.py

@@ -9,10 +9,15 @@
 
 import astropy.units as u
 import numpy as np
+from astropy.coordinates import SkyCoord
+from astropy.io import fits
 from astropy.table import QTable, Table
 from astropy.utils import lazyproperty
+from astropy.wcs import WCS
+from crds import getreferences
 from roman_datamodels.datamodels import ImageModel, MosaicModel
 from roman_datamodels.dqflags import pixel
+from scipy.ndimage import map_coordinates
 
 from romancal import __version__ as romancal_version
 from romancal.skycell import skymap
@@ -138,6 +143,9 @@ class RomanSourceCatalog:
     must also have the same shape and units as the science data array.
     """
 
+    north_galactic_pole_id = "ngp"
+    south_galactic_pole_id = "sgp"
+
     def __init__(
         self,
         model,
@@ -425,6 +433,149 @@ def image_flags(self):
         """
         return np.zeros(self.n_sources, dtype=np.int32)
 
+    @lazyproperty
+    def dust_ebv(self):
+        """
+        Return per-source dust E(B-V) from SFD maps referenced by CRDS.
+
+        This property queries CRDS for north/south SFD map references,
+        evaluates E(B-V) at each source `ra`/`dec` position, and returns
+        one value per source.
+
+        Returns
+        -------
+        result : `~numpy.ndarray`
+            Array of dtype ``float32`` and shape ``(n_sources,)``.
+            If any step fails (CRDS lookup, file I/O, WCS transform, or
+            interpolation), returns all-NaN values with the same shape.
+
+        Warns
+        -----
+        Emits a warning before returning all-NaN values on failure.
+        """
+        try:
+            dust_map_paths = self._get_sfd_map_paths()
+            return self._get_dust_ebv(self.ra, self.dec, dust_map_paths)
+        except Exception as exc:
+            log.warning(
+                f"Failed to compute dust_ebv from SFD reference maps: {exc}; "
+                "setting dust_ebv to NaN."
+            )
+            return np.full(self.n_sources, np.nan, dtype=np.float32)
+
+    def _get_sfd_map_paths(self):
+        """
+        Resolve CRDS references for north/south SFD dust maps.
+
+        Returns
+        -------
+        result : dict
+            Mapping from pole ID to resolved dust-map file path.
+
+        Raises
+        ------
+        Exception
+            Propagates exceptions from CRDS lookups or malformed CRDS return
+            values. The caller (`dust_ebv`) handles these errors and falls
+            back to NaNs.
+        """
+        map_paths = {}
+        for field, pole in (
+            ("north", self.north_galactic_pole_id),
+            ("south", self.south_galactic_pole_id),
+        ):
+            dm_headers = {"roman.meta.instrument.name": "wfi", "roman.field": field}
+            reference = getreferences(
+                dm_headers, reftypes=["dustmap"], observatory="roman"
+            )
+            map_paths[pole] = reference["dustmap"]
+        return map_paths
+
+    @staticmethod
+    def _as_degree_quantity(values):
+        """
+        Convert coordinate-like input to a degree `~astropy.units.Quantity`.
+
+        Parameters
+        ----------
+        values : array-like or `~astropy.units.Quantity`
+            Coordinate values. If already a quantity, converted to degrees.
+
+        Returns
+        -------
+        result : `~astropy.units.Quantity`
+            Values expressed in degree units.
+        """
+        if isinstance(values, u.Quantity):
+            return values.to(u.deg)
+        return np.asarray(values) * u.deg
+
+    def _get_dust_ebv(self, ra, dec, map_paths, order=3):
+        """
+        Interpolate SFD E(B-V) at ICRS coordinates.
+
+        Parameters
+        ----------
+        ra, dec : array-like or `~astropy.units.Quantity`
+            Right ascension and declination. If unitless, interpreted as
+            degrees.
+
+        map_paths : dict
+            Mapping containing FITS file paths for galactic pole ID.
+
+        order : int, optional
+            Interpolation order passed to `scipy.ndimage.map_coordinates`.
+
+        Returns
+        -------
+        result : `~numpy.ndarray`
+            Array of E(B-V) values with the same shape as input positions.
+
+        Raises
+        ------
+        ValueError
+            If ``ra`` and ``dec`` shapes do not match.
+        KeyError
+            If required pole keys are missing from ``map_paths``.
+        Exception
+            Propagates FITS/WCS/interpolation errors to caller, which is
+            handled by `dust_ebv` fallback logic.
+        """
+        ra = self._as_degree_quantity(ra)
+        dec = self._as_degree_quantity(dec)
+
+        ra = np.atleast_1d(ra)
+        dec = np.atleast_1d(dec)
+        if ra.shape != dec.shape:
+            raise ValueError("ra.shape must equal dec.shape")
+
+        if ra.size == 0:
+            return np.array([], dtype=np.float32)
+
+        sky = SkyCoord(ra=ra, dec=dec, frame="icrs")
+        gal = sky.galactic
+        l = gal.l.deg
+        b = gal.b.deg
+
+        out = np.full(ra.shape, np.nan, dtype=np.float32)
+        masks = {
+            self.north_galactic_pole_id: b >= 0,
+            self.south_galactic_pole_id: b < 0,
+        }
+
+        for pole, pole_mask in masks.items():
+            if not np.any(pole_mask):
+                continue
+
+            with fits.open(map_paths[pole]) as hdulist:
+                imwcs = WCS(hdulist[0].header)
+                x, y = imwcs.wcs_world2pix(l[pole_mask], b[pole_mask], 0)
+                out[pole_mask] = map_coordinates(
+                    hdulist[0].data, [y, x], order=order, mode="nearest"
+                ).astype(np.float32)
+
+        return out
+
     def _validate_and_convert_units(self):
         """
         Validate that model data arrays have compatible units and

romancal/source_catalog/tests/test_source_catalog.py

@@ -336,6 +336,25 @@ def test_background(mosaic_model, function_jail):
     assert isinstance(cat, Table)
 
 
+@pytest.mark.parametrize("model_fixture", ("image_model", "mosaic_model"))
+def test_source_catalog_populates_dust_ebv(model_fixture, request):
+    """Ensure prompt source catalogs include a per-source dust_ebv column."""
+    model = request.getfixturevalue(model_fixture)
+    step = SourceCatalogStep()
+    result = step.call(
+        model,
+        bkg_boxsize=50,
+        kernel_fwhm=2.0,
+        snr_threshold=3,
+        npixels=10,
+        save_results=False,
+    )
+    cat = result.source_catalog
+    assert "dust_ebv" in cat.colnames
+    assert len(cat["dust_ebv"]) == len(cat)
+    assert cat["dust_ebv"].dtype == np.float32
+
+
 def test_l2_input_model_unchanged(image_model, function_jail):
     """
     Test that the input model data and error arrays are unchanged after
@@ -572,6 +591,45 @@ def test_l2_source_catalog_keywords(
             assert isinstance(rdm.open(filepath), expected_outputs.get(suffix))
 
 
+@pytest.mark.parametrize(
+    "ra, dec",
+    [
+        (np.array([10.0, 20.0]), np.array([30.0])),
+        (np.array([[10.0, 20.0]]), np.array([30.0, 40.0])),
+    ],
+)
+def test_get_dust_ebv_shape_mismatch_raises(ra, dec):
+    """Raise when RA/Dec input shapes are inconsistent."""
+    cat = object.__new__(RomanSourceCatalog)
+    map_paths = {
+        RomanSourceCatalog.north_galactic_pole_id: "north.fits",
+        RomanSourceCatalog.south_galactic_pole_id: "south.fits",
+    }
+    with pytest.raises(ValueError, match=r"ra\.shape must equal dec\.shape"):
+        cat._get_dust_ebv(ra, dec, map_paths)
+
+
+def test_dust_ebv_property_returns_nan_on_failure(monkeypatch):
+    """Return NaNs when CRDS lookup/interpolation fails."""
+
+    def fail_getreferences(*args, **kwargs):
+        raise RuntimeError()
+
+    monkeypatch.setattr(
+        "romancal.source_catalog._source_catalog.getreferences", fail_getreferences
+    )
+
+    cat = object.__new__(RomanSourceCatalog)
+    cat.ra = np.array([1.0, 2.0, 3.0], dtype=float)
+    cat.dec = np.array([4.0, 5.0, 6.0], dtype=float)
+    cat.n_sources = 3
+
+    result = cat.dust_ebv
+    assert result.dtype == np.float32
+    assert result.shape == (3,)
+    assert np.all(np.isnan(result))
+
+
 @pytest.mark.parametrize(
     "snr_threshold, npixels, nsources, save_results, return_updated_model, expected_result, expected_outputs",
     (