DM-54912: add photometric_scaling to VisitImage

doc/changes/DM-54912.feature.md

@@ -0,0 +1 @@
+Added a `photometric_scaling` component to `VisitImage`.

python/lsst/images/_image.py

@@ -117,7 +117,7 @@ def __init__(
         super().__init__(metadata)
         if isinstance(array_or_fill, np.ndarray):
             if dtype is not None:
-                array = np.array(array_or_fill, dtype=dtype)
+                array = np.array(array_or_fill, dtype=dtype, copy=None)
             else:
                 array = array_or_fill
             if bbox is None:
@@ -483,6 +483,12 @@ def _read_legacy_hdu(
         unit: astropy.units.UnitBase | None = None
         if (fits_unit := hdu.header.pop("BUNIT", None)) is not None:
             unit = astropy.units.Unit(fits_unit, format="fits")
+            if opaque_metadata.get_instrumental_unit() == astropy.units.electron:
+                # Fix incorrect BUNIT='adu' in LSST preliminary_visit_image.
+                if unit == astropy.units.adu:
+                    unit = astropy.units.electron
+                if unit == astropy.units.adu**2:
+                    unit = astropy.units.electron**2
         dx: int = hdu.header.pop("LTV1")
         dy: int = hdu.header.pop("LTV2")
         start = YX(y=-dy, x=-dx)

python/lsst/images/_mask.py

@@ -360,7 +360,7 @@ def __init__(
         if start is not None:
             start = tuple(start)
         if isinstance(array_or_fill, np.ndarray):
-            array = np.array(array_or_fill, dtype=schema.dtype)
+            array = np.array(array_or_fill, dtype=schema.dtype, copy=None)
             if array.ndim != 3:
                 raise ValueError("Mask array must be 3-d.")
             if bbox is None:

python/lsst/images/_masked_image.py

@@ -290,7 +290,7 @@ def read_fits(cls, url: ResourcePathExpression, *, bbox: Box | None = None) -> M
     def from_legacy(
         legacy: Any,
         *,
-        unit: astropy.units.Unit | None = None,
+        unit: astropy.units.UnitBase | None = None,
         plane_map: Mapping[str, MaskPlane] | None = None,
     ) -> MaskedImage:
         """Convert from an `lsst.afw.image.MaskedImage` instance.
@@ -426,13 +426,15 @@ def _read_legacy_hdus(
         hdu_list: astropy.io.fits.HDUList,
         uri: ResourcePathExpression,
         *,
+        opaque_metadata: fits.FitsOpaqueMetadata | None = None,
         preserve_quantization: bool = False,
         plane_map: Mapping[str, MaskPlane] | None = None,
         component: Literal["image", "mask", "variance"] | None,
         fits_wcs_frame: Frame | None = None,
     ) -> Any:
-        opaque_metadata = fits.FitsOpaqueMetadata()
-        opaque_metadata.extract_legacy_primary_header(hdu_list[0].header)
+        if opaque_metadata is None:
+            opaque_metadata = fits.FitsOpaqueMetadata()
+            opaque_metadata.extract_legacy_primary_header(hdu_list[0].header)
         image_bintable_hdu: astropy.io.fits.BinTableHDU | None = None
         variance_bintable_hdu: astropy.io.fits.BinTableHDU | None = None
         result: Any

python/lsst/images/fields/_chebyshev.py

@@ -250,12 +250,32 @@ def _get_archive_tree_type(
 
     @staticmethod
     def from_legacy(
-        legacy: LegacyChebyshevBoundedField, unit: astropy.units.UnitBase | None = None
+        legacy: LegacyChebyshevBoundedField,
+        unit: astropy.units.UnitBase | None = None,
+        bounds: Bounds | None = None,
     ) -> ChebyshevField:
-        """Convert from a legacy `lsst.afw.math.ChebyshevBoundedField`."""
-        return ChebyshevField(
-            bounds=Box.from_legacy(legacy.getBBox()), coefficients=legacy.getCoefficients(), unit=unit
-        )
+        """Convert from a legacy `lsst.afw.math.ChebyshevBoundedField`.
+
+        Parameters
+        ----------
+        legacy
+            Legacy field to convert.
+        unit
+            The units of the returned field (`lsst.afw.math.BoundedField`
+            objects do not know their units).
+        bounds
+            The bounds of the returned field, if they should be different from
+            the bounding box of ``legacy``.
+        """
+        bbox = Box.from_legacy(legacy.getBBox())
+        if bounds is not None:
+            if bounds.bbox != bbox:
+                raise ValueError(
+                    "Custom bounds when converting a ChebyshevBoundedField must not change the bbox."
+                )
+        else:
+            bounds = bbox
+        return ChebyshevField(bounds=bounds, coefficients=legacy.getCoefficients(), unit=unit)
 
     def to_legacy(self) -> LegacyChebyshevBoundedField:
         """Convert to a legacy `lsst.afw.math.ChebyshevBoundedField`."""
@@ -266,9 +286,22 @@ def to_legacy(self) -> LegacyChebyshevBoundedField:
     @staticmethod
     def from_legacy_background(
         legacy_background: LegacyBackground,
+        bounds: Bounds | None = None,
         unit: astropy.units.UnitBase | None = None,
     ) -> ChebyshevField:
-        """Convert from a legacy `lsst.afw.math.BackgroundMI` instance."""
+        """Convert from a legacy `lsst.afw.math.BackgroundMI` instance.
+
+        Parameters
+        ----------
+        legacy
+            Legacy background object to convert.
+        bounds
+            The bounds of the returned field, if they should be different from
+            the bounding box of ``legacy_background``.
+        unit
+            The units of the returned field (`lsst.afw.math.Background`
+            objects do not know their units).
+        """
         from lsst.afw.math import ApproximateControl
 
         approx_control = legacy_background.getBackgroundControl().getApproximateControl()
@@ -281,8 +314,16 @@ def from_legacy_background(
             weight = None
         x = legacy_background.getBinCentersX()
         y = legacy_background.getBinCentersY()
+        bbox = Box.from_legacy(legacy_background.getImageBBox())
+        if bounds is not None:
+            if bounds.bbox != bbox:
+                raise ValueError(
+                    "Custom bounds when converting a Chebyshev background must not change the bbox."
+                )
+        else:
+            bounds = bbox
         return ChebyshevField.fit(
-            Box.from_legacy(legacy_background.getImageBBox()),
+            bounds,
             stats_image.image.array,
             x=x,
             y=y,

python/lsst/images/fields/_concrete.py

@@ -16,27 +16,32 @@
     "FieldSerializationModel",
     "field_from_legacy",
     "field_from_legacy_background",
+    "field_from_legacy_photo_calib",
 )
 
 from typing import TYPE_CHECKING, Annotated, Any
 
 import astropy.units
+import numpy as np
 import pydantic
 
+from .._geom import Bounds
 from ._chebyshev import ChebyshevField, ChebyshevFieldSerializationModel
 from ._product import ProductField, ProductFieldSerializationModel
 from ._spline import SplineField, SplineFieldSerializationModel
 from ._sum import SumField, SumFieldSerializationModel
 
 if TYPE_CHECKING:
     try:
+        from lsst.afw.image import PhotoCalib as LegacyPhotoCalib
         from lsst.afw.math import BackgroundList as LegacyBackgroundList
         from lsst.afw.math import BackgroundMI as LegacyBackground
         from lsst.afw.math import BoundedField as LegacyBoundedField
     except ImportError:
         type LegacyBoundedField = Any  # type: ignore[no-redef]
         type LegacyBackground = Any  # type: ignore[no-redef]
         type LegacyBackgroundList = Any  # type: ignore[no-redef]
+        type LegacyPhotoCalib = Any  # type: ignore[no-redef]
 
 
 # Since Sphinx can't handle doc links to type aliases, whenever we annotate
@@ -61,29 +66,55 @@
 
 
 def field_from_legacy(
-    legacy_bounded_field: LegacyBoundedField, unit: astropy.units.UnitBase | None = None
+    legacy_bounded_field: LegacyBoundedField,
+    bounds: Bounds | None = None,
+    unit: astropy.units.UnitBase | None = None,
 ) -> Field:
     """Convert a legacy `lsst.afw.math.BoundedField` subclass to a `BaseField`
     object.
+
+    Parameters
+    ----------
+    legacy_bounded_field
+        Legacy field to convert.
+    bounds
+        The bounds of the returned field, if they should be different from
+        the bounding box of ``legacy``.
+    unit
+        The units of the returned field (`lsst.afw.math.BoundedField`
+        objects do not know their units).
     """
     from lsst.afw.math import ChebyshevBoundedField, ProductBoundedField
 
     match legacy_bounded_field:
         case ChebyshevBoundedField():
-            return ChebyshevField.from_legacy(legacy_bounded_field, unit=unit)
+            return ChebyshevField.from_legacy(legacy_bounded_field, unit=unit, bounds=bounds)
         case ProductBoundedField():
-            return ProductField.from_legacy(legacy_bounded_field, unit=unit)
+            return ProductField.from_legacy(legacy_bounded_field, unit=unit, bounds=bounds)
         case _:
             raise NotImplementedError(
                 f"Conversion from {type(legacy_bounded_field).__name__} is not supported."
             )
 
 
 def field_from_legacy_background(
-    legacy_background: LegacyBackground | LegacyBackgroundList, unit: astropy.units.UnitBase | None = None
+    legacy_background: LegacyBackground | LegacyBackgroundList,
+    bounds: Bounds | None = None,
+    unit: astropy.units.UnitBase | None = None,
 ) -> Field:
     """Convert a legacy `lsst.afw.math.Background` or
     `lsst.afw.math.BackgroundList` instance to a `BaseField` object.
+
+    Parameters
+    ----------
+    legacy_background
+        Legacy background object to convert.
+    bounds
+        The bounds of the returned field, if they should be different from
+        the bounding box of ``legacy_background``.
+    unit
+        The units of the returned field (`lsst.afw.math.Background`
+        objects do not know their units).
     """
     from lsst.afw.math import ApproximateControl, BackgroundList
 
@@ -95,3 +126,50 @@ def field_from_legacy_background(
         return SplineField.from_legacy_background(legacy_background, unit=unit)
     else:
         return ChebyshevField.from_legacy_background(legacy_background, unit=unit)
+
+
+def field_from_legacy_photo_calib(
+    legacy_photo_calib: LegacyPhotoCalib,
+    bounds: Bounds,
+    instrumental_unit: astropy.units.UnitBase = astropy.units.electron,
+) -> Field | None:
+    """Convert a legacy `lsst.afw.image.PhotoCalib` into a `BaseField` object.
+
+    Parameters
+    ----------
+    legacy_photo_calib
+        Calibration object to convert.
+    bounds
+        Bounds of the returned field.
+    instrumental_unit
+        The instrumental units the legacy calibration transforms from.  These
+        will be used as the denominator of the units of the returned field,
+        with ``astropy.units.nJy`` as the numerator.
+
+    Returns
+    -------
+    `BaseField` | `None`
+        A field that transforms instrumental units to ``nJy``, or `None` if
+        the given calibration object was an identity mapping for a legacy
+        image that already had ``nJy`` pixels.
+    """
+    calibration_mean = legacy_photo_calib.getCalibrationMean()
+    if legacy_photo_calib._isConstant:
+        if calibration_mean == 1.0:
+            # This image's pixels have been calibrated to nJy
+            # already, which means the calibration *from* post-ISR
+            # electrons that we want is stored elsewhere.
+            return None
+        else:
+            return ChebyshevField(
+                bounds,
+                np.array([[calibration_mean]], dtype=np.float64),
+                unit=astropy.units.nJy / instrumental_unit,
+            )
+    else:
+        normalized_field = field_from_legacy(
+            legacy_photo_calib.computeScaledCalibration(),
+            unit=astropy.units.nJy / instrumental_unit,
+            bounds=bounds,
+        )
+        return normalized_field * calibration_mean

python/lsst/images/fields/_product.py

@@ -119,14 +119,28 @@ def _get_archive_tree_type(
 
     @staticmethod
     def from_legacy(
-        legacy: LegacyProductBoundedField, unit: astropy.units.UnitBase | None = None
+        legacy: LegacyProductBoundedField,
+        unit: astropy.units.UnitBase | None = None,
+        bounds: Bounds | None = None,
     ) -> ProductField:
-        """Convert from a legacy `lsst.afw.math.ProductBoundedField`."""
+        """Convert from a legacy `lsst.afw.math.ProductBoundedField`.
+
+        Parameters
+        ----------
+        legacy
+            Legacy field to convert.
+        unit
+            The units of the returned field (`lsst.afw.math.BoundedField`
+            objects do not know their units).
+        bounds
+            The bounds of the returned field, if they should be different from
+            the bounding box of ``legacy``.
+        """
         from ._concrete import field_from_legacy
 
         legacy_factors = legacy.getFactors()
-        operands = [field_from_legacy(f) for f in legacy_factors[:-1]]
-        operands.append(field_from_legacy(legacy_factors[-1], unit=unit))
+        operands = [field_from_legacy(f, bounds=bounds) for f in legacy_factors[:-1]]
+        operands.append(field_from_legacy(legacy_factors[-1], unit=unit, bounds=bounds))
         return ProductField(operands)
 
     def to_legacy(self) -> LegacyProductBoundedField:

python/lsst/images/fields/_spline.py

@@ -191,10 +191,22 @@ def _get_archive_tree_type(
     @staticmethod
     def from_legacy_background(
         legacy_background: LegacyBackground,
+        bounds: Bounds | None = None,
         unit: astropy.units.UnitBase | None = None,
     ) -> SplineField:
         """Convert from a legacy `lsst.afw.math.BackgroundMI` instance.
 
+        Parameters
+        ----------
+        legacy
+            Legacy background object to convert.
+        bounds
+            The bounds of the returned field, if they should be different from
+            the bounding box of ``legacy_background``.
+        unit
+            The units of the returned field (`lsst.afw.math.Background`
+            objects do not know their units).
+
         Notes
         -----
         `SplineField.render` and the `lsst.afw` background interpolator both
@@ -216,7 +228,11 @@ def from_legacy_background(
         x = legacy_background.getBinCentersX()
         y = legacy_background.getBinCentersY()
         return SplineField(
-            Box.from_legacy(legacy_background.getImageBBox()), stats_image.image.array, x=x, y=y, unit=unit
+            Box.from_legacy(legacy_background.getImageBBox()) if bounds is None else bounds,
+            stats_image.image.array,
+            x=x,
+            y=y,
+            unit=unit,
         )
 
     def _make_1d_interpolator(self, loc: np.ndarray, val: np.ndarray) -> Akima1DInterpolator | None:

python/lsst/images/fields/_sum.py

@@ -126,12 +126,27 @@ def _get_archive_tree_type(
     @staticmethod
     def from_legacy_background(
         legacy_background: LegacyBackgroundList,
+        bounds: Bounds | None = None,
         unit: astropy.units.UnitBase | None = None,
     ) -> SumField:
-        """Convert from a legacy `lsst.afw.math.BackgroundList` instance."""
+        """Convert from a legacy `lsst.afw.math.BackgroundList` instance.
+
+        Parameters
+        ----------
+        legacy
+            Legacy background object to convert.
+        bounds
+            The bounds of the returned field, if they should be different from
+            the bounding box of ``legacy_background``.
+        unit
+            The units of the returned field (`lsst.afw.math.BackgroundList`
+            objects do not know their units).
+        """
         from ._concrete import field_from_legacy_background
 
-        return SumField([field_from_legacy_background(b, unit) for b, *_ in legacy_background])
+        return SumField(
+            [field_from_legacy_background(b, bounds=bounds, unit=unit) for b, *_ in legacy_background]
+        )
 
 
 class SumFieldSerializationModel(ArchiveTree):

python/lsst/images/fits/_common.py

@@ -383,6 +383,15 @@ def subset(self, bbox: Box) -> FitsOpaqueMetadata:
         # Docstring inherited.
         return FitsOpaqueMetadata(headers=self.headers)
 
+    def get_instrumental_unit(self) -> astropy.units.UnitBase | None:
+        """Extract the ``LSST ISR UNIT`` key from the primary header (if it
+        exists) and wrap it with Astropy.
+        """
+        if (primary_header := self.headers.get(ExtensionKey())) is not None:
+            if (instrumental_unit_str := primary_header.get("LSST ISR UNITS")) is not None:
+                return astropy.units.Unit(instrumental_unit_str)
+        return None
+
 
 def add_offset_wcs(header: astropy.io.fits.Header, *, x: int | float, y: int | float, key: str = "A") -> None:
     """Add a trivial FITS WCS to a header that applies the appropriate offset