Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Restore flux units to source catalog table #1512

Merged
merged 3 commits into from
Nov 13, 2024
Merged

Restore flux units to source catalog table #1512

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
2bc68ed
Restore flux units to source catalog table
larrybradley Nov 13, 2024
ce47849
Add changelog entry
larrybradley Nov 13, 2024
109d7eb
Merge branch 'main' into sourcecat-units
schlafly Nov 13, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions changes/1512.source_catalog.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
Restored flux units in source catalog table.
63 changes: 34 additions & 29 deletions romancal/source_catalog/source_catalog.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -98,13 +98,8 @@ def __init__(
self.convolved_data = convolved_data
self.aperture_params = aperture_params
self.kernel_sigma = kernel_fwhm * gaussian_fwhm_to_sigma
self.flux_unit = flux_unit
self.fit_psf = fit_psf

self.sb_unit = "MJy/sr"
self.l2_unit = "DN/s"
self.l2_conv_factor = self.model.meta.photometry.conversion_megajanskys

if len(ci_star_thresholds) != 2:
raise ValueError("ci_star_thresholds must contain only 2 items")
self.ci_star_thresholds = ci_star_thresholds
Expand All @@ -122,6 +117,12 @@ def __init__(
self.get_psf_photometry_catalog_colnames_mapping()
)

self.flux_unit = flux_unit
self.l2_conv_factor = self.model.meta.photometry.conversion_megajanskys
self.sb_to_flux = (1.0 * (u.MJy / u.sr) * self.pixel_area).to(
u.Unit(self.flux_unit)
)

@lazyproperty
def _pixscale_angle(self):
ysize, xsize = self.model.data.shape
Expand Down Expand Up @@ -154,20 +155,20 @@ def pixel_area(self):
"""
The pixel area in steradians.

If the value is a negative placeholder (e.g., -999999 sr), the
value is calculated from the WCS at the center of the image.
If the meta.photometry.pixel_area value is a negative
placeholder (e.g., -999999 sr), the value is calculated from the
WCS at the center of the image.
"""
pixel_area = self.model.meta.photometry.pixel_area
pixel_area = self.model.meta.photometry.pixel_area * u.sr
if pixel_area < 0:
pixel_area = (self._pixel_scale**2).to(u.sr).value
pixel_area = (self._pixel_scale**2).to(u.sr)
return pixel_area

def convert_l2_to_sb(self):
"""
Convert a level-2 image from units of DN/s to MJy/sr (surface
brightness).
"""

# the conversion in done in-place to avoid making copies of the data;
# use a dictionary to set the value to avoid on-the-fly validation
self.model["data"] *= self.l2_conv_factor
Expand All @@ -181,7 +182,6 @@ def convert_sb_to_l2(self):

This is the inverse operation of `convert_l2_to_sb`.
"""

# the conversion in done in-place to avoid making copies of the data;
# use a dictionary to set the value to avoid on-the-fly validation
self.model["data"] /= self.l2_conv_factor
Expand All @@ -193,28 +193,34 @@ def convert_sb_to_flux_density(self):
Convert the data and error arrays from MJy/sr (surface
brightness) to flux density units.

The arrays are converted in-place to Quantity arrays.

The flux density unit is defined by self.flux_unit.
"""

# the conversion in done in-place to avoid making copies of the data;
# use a dictionary to set the value to avoid on-the-fly validation
self.model["data"] *= self.pixel_area

self.model["err"] *= self.pixel_area

self.convolved_data *= self.pixel_area
self.model["data"] *= self.sb_to_flux.value
self.model["err"] *= self.sb_to_flux.value
self.convolved_data *= self.sb_to_flux.value
self.model["data"] <<= self.sb_to_flux.unit
self.model["err"] <<= self.sb_to_flux.unit
self.convolved_data <<= self.sb_to_flux.unit

def convert_flux_density_to_sb(self):
"""
Convert the data and error arrays from flux density units to
MJy/sr (surface brightness).

The units are also removed from the arrays.

This is the inverse operation of `convert_sb_to_flux_density`.
"""

self.model["data"] /= self.pixel_area
self.model["err"] /= self.pixel_area
self.convolved_data /= self.pixel_area
self.model["data"] /= self.sb_to_flux
self.model["err"] /= self.sb_to_flux
self.convolved_data /= self.sb_to_flux
self.model["data"] = self.model["data"].value
self.model["err"] = self.model["err"].value
self.convolved_data = self.convolved_data.value

def convert_flux_to_abmag(self, flux, flux_err):
"""
Expand All @@ -234,15 +240,11 @@ def convert_flux_to_abmag(self, flux, flux_err):
with warnings.catch_warnings():
warnings.simplefilter("ignore", category=RuntimeWarning)

# exact AB mag zero point
flux_zpt = 10 ** (-0.4 * 48.60)

abmag_zpt = 2.5 * np.log10(flux_zpt)
abmag = -2.5 * np.log10(flux) + abmag_zpt
abmag = flux.to(u.ABmag).value
abmag_err = 2.5 * np.log10(1.0 + (flux_err / flux))

# handle negative fluxes
idx = flux < 0
idx = flux.value < 0
abmag[idx] = np.nan
abmag_err[idx] = np.nan

Expand Down Expand Up @@ -534,7 +536,7 @@ def _aper_local_background(self):
bkg_median = []
bkg_std = []
for mask in bkg_aper_masks:
bkg_data = mask.get_values(self.model.data)
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))
Expand All @@ -545,6 +547,9 @@ def _aper_local_background(self):
# 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
Expand Down Expand Up @@ -754,7 +759,7 @@ def _daofind_convolved_data(self):
The DAOFind convolved data.
"""
return ndimage.convolve(
self.model.data, self._daofind_kernel, mode="constant", cval=0.0
self.model.data.value, self._daofind_kernel, mode="constant", cval=0.0
)

@lazyproperty
Expand Down