JP-4161: Fix SourceCatalogStep fluxes failing to convert from Mjy/sr to Jy (#9939)

Author: emolter

changes/9939.source_catalog.rst

@@ -0,0 +1,2 @@
+Fix a bug where fluxes were not being converted from Mjy/sr to Jy.
+Fix a bug where fluxes included background flux when they shouldn't.

jwst/source_catalog/source_catalog_step.py

@@ -104,6 +104,9 @@ def process(self, input_model):
 
             coverage_mask = np.isnan(model.err) | (model.wht == 0)
 
+            # convert to Jy before calling make_tweakreg_catalog so the outputs end up in Jy
+            JWSTSourceCatalog.convert_mjysr_to_jy(model)
+
             starfinder_kwargs = {
                 "sigma_radius": self.sigma_radius,
                 "minsep_fwhm": self.minsep_fwhm,
@@ -139,7 +142,6 @@ def process(self, input_model):
                 log.warning("No sources found in the image. Catalog will be empty.")
                 return None
 
-            JWSTSourceCatalog.convert_mjysr_to_jy(model)
             ci_star_thresholds = (self.ci1_star_threshold, self.ci2_star_threshold)
             catobj = JWSTSourceCatalog(
                 model,

jwst/source_catalog/tests/test_source_catalog.py

@@ -3,6 +3,7 @@
 import numpy as np
 import pytest
 import stdatamodels.jwst.datamodels as dm
+from astropy.modeling import models
 from astropy.table import QTable
 from numpy.testing import assert_allclose
 from photutils.datasets import make_gwcs
@@ -29,7 +30,7 @@ def nircam_model():
     model = dm.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.photometry.pixelarea_steradians = 1.0e-13
     model.meta.wcs = make_gwcs(data.shape)
     model.meta.wcsinfo = {
         "ctype1": "RA---TAN",
@@ -80,7 +81,7 @@ def nircam_model_without_apcorr():
     model = dm.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.photometry.pixelarea_steradians = 1.0e-13
     model.meta.wcs = make_gwcs(data.shape)
     model.meta.wcsinfo = {
         "ctype1": "RA---TAN",
@@ -130,20 +131,22 @@ def test_source_catalog(nircam_model, npixels, nsources):
 
     if npixels == 5 and nsources == 2:
         # test values of some specific computed quantities
-        assert np.isclose(cat["xcentroid"][1], 19.46399720865899)
-        assert np.isclose(cat["ycentroid"][1], 41.95288393407728)
-        assert np.isclose(cat["aper_bkg_flux"][1].value, 1400000.0)
-        assert np.isclose(cat["aper_bkg_flux_err"][1].value, 85223.70700074881)
-        assert np.isclose(cat["CI_50_30"][1], 2.3342599432074653)
+        assert np.isclose(cat["xcentroid"][1], 19.453064764431833)
+        assert np.isclose(cat["ycentroid"][1], 41.963065678485115)
+        assert np.isclose(cat["aper_bkg_flux"][1].value, 1.40e-7)
+        assert np.isclose(cat["aper_bkg_flux_err"][1].value, 8.52237054e-09)
+        assert np.isclose(cat["CI_50_30"][1], 2.352272196434725)
         assert np.isclose(cat["sharpness"][1], 0.9102634628764403)
         assert np.isclose(cat["roundness"][1], 1.5954264)
-        assert np.isclose(cat["nn_dist"][1].value, 53.0737632103816)
-        assert np.isclose(cat["isophotal_flux"][1], 930.9999841451645)
-        assert np.isclose(cat["isophotal_flux_err"][1], 3.6102633)
-        assert np.isclose(cat["semimajor_sigma"][1].value, 18.847635525516534)
-        assert np.isclose(cat["semiminor_sigma"][1].value, 7.031371175038476)
-        assert np.isclose(cat["ellipticity"][1], 0.626936165784871)
-        assert np.isclose(cat["orientation"][1].value, -72.75413766990114)
+        assert np.isclose(cat["nn_dist"][1].value, 53.07168773319497)
+        assert np.isclose(cat["isophotal_flux"][1].value, 7.940753383195442e-05)
+        assert cat["isophotal_flux_err"][1].unit == "Jy"
+        assert np.isclose(cat["isophotal_flux_err"][1].value, 3.6102634e-07)
+        assert np.isclose(cat["isophotal_abmag"][1], 19.150345729246215)
+        assert np.isclose(cat["semimajor_sigma"][1].value, 18.84372169911978)
+        assert np.isclose(cat["semiminor_sigma"][1].value, 7.024931388267103)
+        assert np.isclose(cat["ellipticity"][1].value, 0.62720043)
+        assert np.isclose(cat["orientation"][1].value, -72.78329207140818)
 
 
 def test_source_catalog_no_sources(nircam_model, monkeypatch):
@@ -191,14 +194,20 @@ def test_source_catalog_point_sources(finder, nircam_model, tmp_cwd):
     """Test the three source finding algorithms with point sources."""
     data = np.random.default_rng(seed=123).normal(0, 0.5, size=(101, 101))
 
-    # make a point source with some size that looks a bit like a psf, no need to be realistic
-    point_source = np.ones((7, 7))
-    point_source[1:6, 1:6] = 3.0
-    point_source[2:5, 2:5] = 5.0
-    point_source[3, 3] = 10.0
+    # Create coordinate grids for Airy disk models
+    y, x = np.mgrid[0:101, 0:101]
 
-    data[30:37, 30:37] = point_source
-    data[70:77, 70:77] = point_source
+    # Create two Airy disk point sources
+    # First source at position (33, 33)
+    airy1 = models.AiryDisk2D(amplitude=10.0, x_0=33, y_0=33, radius=5.0)
+    source1 = airy1(x, y)
+
+    # Second source at position (73, 73)
+    airy2 = models.AiryDisk2D(amplitude=10.0, x_0=73, y_0=73, radius=5.0)
+    source2 = airy2(x, y)
+
+    # Add the sources to the background
+    data += source1 + source2
 
     nircam_model.data = data
 
@@ -210,6 +219,8 @@ def test_source_catalog_point_sources(finder, nircam_model, tmp_cwd):
         npixels=5,
         bkg_boxsize=50,
         kernel_fwhm=2.0,
+        roundlo=-1.0,
+        sharplo=0.0,
         starfinder=finder,
         save_results=True,
     )

jwst/tweakreg/tweakreg_catalog.py

@@ -412,7 +412,8 @@ def make_tweakreg_catalog(
         with warnings.catch_warnings():
             # suppress warning about NaNs being automatically masked - this is desired
             warnings.simplefilter("ignore", AstropyUserWarning)
-            threshold_img = bkg.background + (snr_threshold * bkg.background_rms)
+            threshold_img = snr_threshold * bkg.background_rms
+        data = model.data - bkg.background
     except ValueError as e:
         log.warning(f"Error determining sky background: {e.args[0]}")
         sources = _empty_table()
@@ -425,7 +426,7 @@ def make_tweakreg_catalog(
             "ignore", category=NoDetectionsWarning, message="No sources were found"
         )
         sources, segmentation_image = starfinder(
-            model.data,
+            data,
             threshold_img,
             kernel_fwhm,
             mask=coverage_mask,