[05/29] Preserve arm exposures, charge binning and independent detector readouts

scopesim/detector/detector.py

@@ -25,6 +25,7 @@
 
 class Detector:
     def __init__(self, header, cmds=None, **kwargs):
+        self._native_header = header.copy()
         self._hdu = ImageHDU(header=header, data=zeros_from_header(header))
         self.meta = {}
         self.meta.update(header)
@@ -51,10 +52,9 @@ def extract_from(self, image_plane, spline_order=1, reset=True):
             self.header.update(sky_wcs.to_header())
 
     def reset(self):
-        """Reset internal HDU data to all-zeros-array."""
-        # The detector might have been converted to integers by the
-        # Quantization effect, so it is not possible to use zeros_like.
-        self._hdu.data = np.zeros(self._hdu.data.shape)
+        """Reset charge and restore the configured, unbinned detector geometry."""
+        self._hdu = ImageHDU(header=self._native_header, data=zeros_from_header(self._native_header))
+        self.meta.update(self._native_header)
 
     @property
     def hdu(self):

scopesim/effects/detector_list.py

@@ -330,16 +330,17 @@ def detector_headers(self, ids=None):
                 hdr["PC2_1"] = -sang
                 hdr["PC2_2"] = cang
 
-            # hdr["GAIN"] = row["gain"]
-            if "id" in row:
-                hdr["DET_ID"] = row["id"]
-                hdr["EXTNAME"] = f"DET_{row['id']}"
-
             row_dict = {
                 f"HIERARCH {col.upper()}" if len(col) > 8 else col: row[col]
                 for col in row.colnames
             }
             hdr.update(row_dict)
+            if "id" in row.colnames:
+                # Keep routing identifiers numeric when detector metadata is
+                # serialized to FITS (NumPy integers otherwise become strings).
+                hdr["ID"] = int(row["id"])
+                hdr["DET_ID"] = int(row["id"])
+                hdr["EXTNAME"] = f"DET_{row['id']}"
             hdrs.append(hdr)
 
         return hdrs

scopesim/effects/electronic/__init__.py

@@ -52,9 +52,9 @@ def apply_to(self, obj, **kwargs):
 
 from .electrons import LinearityCurve, ADConversion, InterPixelCapacitance
 from .noise import (Bias, PoorMansHxRGReadoutNoise, BasicReadoutNoise,
-                    PixelResponseNonUniformity, ShotNoise, DarkCurrent)
+                    PixelResponseNonUniformity, ShotNoise, DarkCurrent, ReadoutNoiseFromTable)
 from .exposure import AutoExposure, ExposureIntegration, ExposureOutput
-from .pixels import ReferencePixelBorder, BinnedImage, UnequalBinnedImage
+from .pixels import ReferencePixelBorder, BinnedImage, UnequalBinnedImage, OnChipBinning
 from .dmps import DetectorModePropertiesSetter
 
 

scopesim/effects/electronic/dmps.py

@@ -23,6 +23,9 @@ class DetectorModePropertiesSetter(Effect):
     mode_properties : dict
         A dictionary containing the DET parameters to be changed for each mode.
         See below for an example yaml entry.
+    mode_key : str
+        Command selecting the mode, by default ``!OBS.detector_readout_mode``.
+        Separate keys allow independent detector arms to select their modes.
 
     Examples
     --------
@@ -70,21 +73,20 @@ def __init__(self, **kwargs):
         check_keys(self.meta, self.required_keys, action="error")
 
         self.mode_properties = kwargs["mode_properties"]
+        self.mode_key = kwargs.get("mode_key", "!OBS.detector_readout_mode")
 
     def apply_to(self, obj, **kwargs):
         logger.debug("Entering DetectorModePropertiesSetter.apply_to")
-        mode_name = kwargs.get(
-            "detector_readout_mode",
-            from_currsys("!OBS.detector_readout_mode", self.cmds),
-        )
-
+        mode_name = kwargs.get("detector_readout_mode",
+                               from_currsys(self.mode_key,
+                                            self.cmds))
         if isinstance(obj, ImagePlane) and mode_name == "auto":
             mode_name = self.select_mode(obj, **kwargs)
             logger.info("Detector mode set to %s", mode_name)
 
         self.meta["detector_readout_mode"] = mode_name
         props_dict = self.mode_properties[mode_name]
-        self.cmds["!OBS.detector_readout_mode"] = mode_name
+        self.cmds[self.mode_key] = mode_name
         for key, value in props_dict.items():
             self.cmds[key] = value
             logger.debug("DMPS setting %s: %s", key, value)

scopesim/effects/electronic/electrons.py

@@ -263,7 +263,7 @@ def __init__(self, **kwargs):
         super().__init__(**kwargs)
         params = {
             "dtype": "uint16",
-            "gain": 1.      # default, usually overridden from yaml
+            "gain": "!DET.gain"
         }
         self.meta.update(params)
         self.meta.update(kwargs)
@@ -297,17 +297,18 @@ def _get_gain(self, det_id: int) -> Real:
         # Apply the gain value (copy from DarkCurrent)
         # Note that this does not cater for the case where the gain is given
         # as a plain dictionary. Should we implement that?
-        if hasattr(self.cmds["!DET.gain"], "dic"):
-            gain = self.cmds["!DET.gain"].dic[det_id]
+        configured_gain = from_currsys(self.meta["gain"], self.cmds)
+        if hasattr(configured_gain, "dic"):
+            gain = configured_gain.dic[det_id]
             logger.info(f"Detector {det_id}: applying gain {gain}")
             return gain
-        if isinstance(self.cmds["!DET.gain"], Real):
-            gain = self.cmds["!DET.gain"]
+        if isinstance(configured_gain, Real):
+            gain = configured_gain
             logger.info(f"Applying gain {gain}")
             return gain
         raise ValueError(
             f"{self.__class__.__name__}.meta['gain'] must be either "
-            f"dict or float, but is {self.cmds['!DET.gain']}")
+            f"dict or float, but is {configured_gain}")
 
     def _apply_to_det(self, det: Detector) -> None:
         """Subclasses can override if more params needed in call."""

scopesim/effects/electronic/exposure.py

@@ -287,6 +287,9 @@ def apply_to(self, obj, **kwargs):
                 "readout call."
             )
 
+        obj.header["EXPTIME"] = (dit * ndit, "Total accumulated integration time [s]")
+        obj.header['DIT'] = (dit, 'Integration time per exposure [s]')
+        obj.header['NDIT'] = ndit
         obj.data = obj.data * dit * ndit
 
         return obj

scopesim/effects/electronic/noise.py

@@ -54,7 +54,7 @@ def __call__(
         return data + dark_level * dit * ndit
 
     def _get_dark_level(self, det_id: int) -> float:
-        dark_level = float(from_currsys(self.meta["value"], self.cmds))
+        dark_level = from_currsys(self.meta["value"], self.cmds)
         if isinstance(dark_level, Real):
             return dark_level
         if isinstance(dark_level, Mapping):
@@ -286,6 +286,23 @@ def _pseudo_random_field(
         return image
 
 
+class ReadoutNoiseFromTable(BasicReadoutNoise):
+    """Interpolate measured read noise in electrons as a function of DIT [s]."""
+
+    required_keys = {"filename", "dit", "ndit"}
+
+    def apply_to(self, det, **kwargs):
+        if not isinstance(det, Detector):
+            return det
+        dit = from_currsys(self.meta["dit"], self.cmds)
+        times = self.table["dit"]
+        if not times[0] <= dit <= times[-1]:
+            raise ValueError(f"DIT {dit} s is outside the measured read-noise table")
+        self.meta["noise_std"] = np.interp(dit, times, self.table["read_noise"])
+        det.header["RDNOISE"] = (self.meta["noise_std"], "Read noise per integration [electron]")
+        return super().apply_to(det, **kwargs)
+
+
 class PixelResponseNonUniformity(RandomEffect):
     """Pixel Response Non-Uniformity (PRNU).
 

scopesim/effects/electronic/pixels.py

@@ -93,19 +93,55 @@ class BinnedImageBase(ElectronicEffect):
 
     z_order: ClassVar[tuple[int, ...]] = (870,)
 
+    def __call__(self, data: ArrayLike) -> NDArray:
+        bx, by = self.bin_factors
+        ny, nx = data.shape
+        return data.reshape(ny//by, by, nx//bx, bx).sum(axis=(1, 3))
+
+    def bin_header(self, native_header):
+        """WCS of the summed pixel blocks, also used by trace coordinates."""
+        bx, by = self.bin_factors
+        if not isinstance(bx, int) or not isinstance(by, int) or min(bx, by) < 1:
+            raise ValueError('Bin factors must be positive integers')
+        if 'allowed_bins' in self.meta and [bx, by] not in self.meta['allowed_bins']:
+            raise ValueError(f'Unsupported binning: {bx} x {by}')
+        header = native_header.copy()
+        for axis, factor in ((1, bx), (2, by)):
+            if header[f'NAXIS{axis}'] % factor:
+                raise ValueError('Binning must divide the full detector dimensions')
+            header[f'NAXIS{axis}'] //= factor
+        for suffix in [''] + list('ABCDEFGHIJKLMNOPQRSTUVWXYZ'):
+            if f'CRPIX1{suffix}' not in header:
+                continue
+            for axis, factor in ((1, bx), (2, by)):
+                header[f'CRPIX{axis}{suffix}'] = (header[f'CRPIX{axis}{suffix}'] - .5)/factor + .5
+            if f'CD1_1{suffix}' in header:
+                for i in (1, 2):
+                    for j, factor in ((1, bx), (2, by)):
+                        header[f'CD{i}_{j}{suffix}'] *= factor
+            else:
+                header[f'CDELT1{suffix}'] *= bx
+                header[f'CDELT2{suffix}'] *= by
+                if f'PC1_2{suffix}' in header:
+                    header[f'PC1_2{suffix}'] *= by/bx
+                if f'PC2_1{suffix}' in header:
+                    header[f'PC2_1{suffix}'] *= bx/by
+        header['BINX'], header['BINY'] = bx, by
+        return header
+
+    def _apply_to_det(self, det):
+        header = self.bin_header(det.header)
+        det.data = self(det.data)
+        det.header.update(header)
+        det.meta.update(header)
+
 
 class BinnedImage(BinnedImageBase):
     required_keys = {"bin_size"}
 
-    def __call__(self, data: ArrayLike) -> NDArray:
-        height, width = data.shape
-        binned_data = data.reshape((
-            height//self.bin_size,
-            self.bin_size,
-            width//self.bin_size,
-            self.bin_size
-        )).sum(axis=3).sum(axis=1)
-        return binned_data
+    @property
+    def bin_factors(self):
+        return self.bin_size, self.bin_size
 
     @property
     def bin_size(self) -> int:
@@ -115,15 +151,9 @@ def bin_size(self) -> int:
 class UnequalBinnedImage(BinnedImageBase):
     required_keys = {"binx","biny"}
 
-    def __call__(self, data: ArrayLike) -> NDArray:
-        height, width = data.shape
-        binned_data = data.reshape((
-            height//self.biny,
-            self.biny,
-            width//self.binx,
-            self.binx
-        )).sum(axis=3).sum(axis=1)
-        return binned_data
+    @property
+    def bin_factors(self):
+        return self.binx, self.biny
 
     @property
     def binx(self) -> int:
@@ -132,3 +162,14 @@ def binx(self) -> int:
     @property
     def biny(self) -> int:
         return from_currsys(self.meta["biny"], self.cmds)
+
+
+class OnChipBinning(UnequalBinnedImage):
+    """Sum charge in FITS x/y pixel blocks before adding read noise.
+
+    Unlike a software sum after readout, each output pixel receives one read
+    noise draw. Configure this effect between charge integration and read noise.
+    WCS reference pixels refer to the centers of the new pixel blocks.
+    """
+
+    z_order = (815,)

scopesim/effects/fits_headers.py

@@ -153,6 +153,9 @@ class ExtraFitsKeywords(Effect):
           header['HIERARCH ESO DET DIT'] = (5, '[s] exposure length')
           header['EXTNAME'] = "DET2.DATA"
 
+    A group may additionally specify ``detector_id``. It then applies only to
+    HDULists containing that ``DET_ID``, including their primary headers.
+
     Resolved and un-resolved keywords
     ---------------------------------
     ScopeSim uses bang-strings to resolve global parameters.
@@ -285,6 +288,12 @@ def apply_to(self, hdul, **kwargs):
             return hdul
 
         for dic in self.dict_list:
+            # Each detector array produces its own HDUList. A detector selector
+            # also permits arm-specific metadata in that array's primary HDU.
+            if "detector_id" in dic:
+                detector_ids = [hdu.header["DET_ID"] for hdu in hdul if "DET_ID" in hdu.header]
+                if dic["detector_id"] not in detector_ids:
+                    continue
             resolved = flatten_dict(dic.get("keywords", {}), resolve=True,
                                     optics_manager=opt_train)
             unresolved = flatten_dict(dic.get("unresolved_keywords", {}))

scopesim/effects/mosaic_trace_list.py

@@ -255,6 +255,9 @@ def apply_to(self, det, **kwargs):
 
         image = det.data
         detwcs = WCS(det.header, key="D")
+        # ExposureIntegration has already recorded the actual accumulated time.
+        # New table HDUs need those cards, but not the detector-image WCS.
+        exposure_header = fits.Header(card for card in det.header.cards if card.keyword in ("DIT", "NDIT", "EXPTIME"))
 
         output_format = from_currsys(self.meta["format"], self.cmds)
 
@@ -309,6 +312,7 @@ def apply_to(self, det, **kwargs):
         else:
             logger.error("Unknown output format: %s", output_format)
 
+        det.header.update(exposure_header)
         return det