Merge pull request #392 from LCOGT/stacked-images

Added support for stacked images

.github/workflows/e2e.yaml

@@ -43,9 +43,21 @@ jobs:
           # Wait for banzai to be ready
           kubectl wait --for=condition=Ready --timeout=60m pod/banzai-e2e-test
 
-      - name: Test Master Bias Creation
+      - name: Test Super Bias Creation
         run: |
-          kubectl exec banzai-e2e-test -c banzai-listener -- pytest -o log_cli=true -s --pyargs banzai --durations=0 --junitxml=/archive/engineering/pytest-master-bias.xml -m master_bias
+          kubectl exec banzai-e2e-test -c banzai-listener -- pytest -s --pyargs banzai --durations=0 --junitxml=/archive/engineering/pytest-master-bias.xml -m master_bias
+      
+      - name: Test Super Dark Creation
+        run: |
+          kubectl exec banzai-e2e-test -c banzai-listener -- pytest -s --pyargs banzai --durations=0 --junitxml=/archive/engineering/pytest-master-dark.xml -m master_dark
+      
+      - name: Test Super Flat Creation
+        run: |
+          kubectl exec banzai-e2e-test -c banzai-listener -- pytest -s --pyargs banzai --durations=0 --junitxml=/archive/engineering/pytest-master-flat.xml -m master_flat
+      
+      - name: Test Science Frame Creation
+        run: |
+          kubectl exec banzai-e2e-test -c banzai-listener -- pytest -s --pyargs banzai --durations=0 --junitxml=/archive/engineering/pytest-science-files.xml -m science_files
       
       - name: Cleanup
         run: |

.github/workflows/unit-tests.yml

@@ -18,11 +18,6 @@ jobs:
             python-version: 3.9
             toxenv: build_docs
 
-          - name: Python 3.9 with developer version of astropy and numpy
-            os: ubuntu-latest
-            python-version: 3.9
-            toxenv: py39-test-devdeps
-
           - name: Python 3.8 with minimal dependencies
             os: ubuntu-latest
             python-version: '3.8'

CHANGES.md

@@ -1,18 +1,22 @@
-1.17.0 (2023-04-24)
+1.18.0 (2024-08-09)
+- Added support for frames that are composed of sub-exposures that 
+  are stacked at site
+
+1.17.0 (2024-04-24)
 -------------------
 - We now omit sources in the photometry stage that have an area larger than 1000 pixels as they lead to long 
   processing times and are almost invariably spurious.
 
-1.16.1 (2023-04-23)
+1.16.1 (2024-04-23)
 -------------------
 - Correction to aperture photometry. We were incorrectly using the radius instead of the diameter
 
-1.16.0 (2023-04-18)
+1.16.0 (2024-04-18)
 -------------------
 - Calibration frames are now associated with output data products rather than frames
   so that we have more than one calibration data product produced per frame.
 
-1.15.2 (2023-04-12)
+1.15.2 (2024-04-12)
 -------------------
 - Fix to fpacking data when the image data array is None
 

Dockerfile

@@ -24,3 +24,5 @@ COPY --chown=10087:10000 . /lco/banzai
 ENV PATH /home/archive/envs/banzai/bin:$PATH
 
 RUN /home/archive/envs/banzai/bin/pip install --no-cache-dir /lco/banzai/ 
+
+RUN cp /lco/banzai/pytest.ini /home/archive/pytest.ini

banzai/bias.py

@@ -17,6 +17,8 @@ def calibration_type(self):
         return 'BIAS'
 
     def make_master_calibration_frame(self, images):
+        for image in images:
+            image /= image.n_sub_exposures
         master_image = super(BiasMaker, self).make_master_calibration_frame(images)
         master_image.bias_level = np.mean([image.bias_level for image in images if image.bias_level is not None])
         return master_image
@@ -31,9 +33,9 @@ def calibration_type(self):
         return 'bias'
 
     def apply_master_calibration(self, image, master_calibration_image):
-        image -= master_calibration_image.bias_level
+        image -= master_calibration_image.bias_level * image.n_sub_exposures
         image.meta['BIASLVL'] = master_calibration_image.bias_level, 'Bias level that was removed after overscan'
-        image -= master_calibration_image
+        image -= master_calibration_image * image.n_sub_exposures
         image.meta['L1IDBIAS'] = master_calibration_image.filename, 'ID of bias frame'
         image.meta['L1STATBI'] = 1, "Status flag for bias frame correction"
         return image
@@ -65,7 +67,9 @@ def __init__(self, runtime_context):
     def do_stage(self, image):
         bias_level = stats.sigma_clipped_mean(image.data, 3.5, mask=image.mask)
         image -= bias_level
-        image.meta['BIASLVL'] = bias_level, 'Bias level that was removed after overscan'
+        # This is only run for bias frames and n_sub_exposures should always be 1
+        # but we divide it here for consistency to cover pathological uses cases
+        image.meta['BIASLVL'] = bias_level / image.n_sub_exposures, 'Bias level that was removed after overscan'
         return image
 
 

banzai/data.py

@@ -132,7 +132,8 @@ def __init__(self, data: Union[np.array, Table], meta: fits.Header,
                  mask: np.array = None, name: str = '', uncertainty: np.array = None, memmap=True):
         super().__init__(data=data, meta=meta, mask=mask, name=name, memmap=memmap)
         if uncertainty is None:
-            uncertainty = self.read_noise * np.ones(data.shape, dtype=data.dtype) / self.gain
+            uncertainty = self.read_noise * np.sqrt(self.n_sub_exposures) * np.ones(data.shape, dtype=data.dtype)
+            uncertainty /= self.gain
         self.uncertainty = self._init_array(uncertainty)
         self._detector_section = Section.parse_region_keyword(self.meta.get('DETSEC'))
         self._data_section = Section.parse_region_keyword(self.meta.get('DATASEC'))
@@ -155,6 +156,15 @@ def __imul__(self, value):
         self.meta['MAXLIN'] *= value
         return self
 
+    def __mul__(self, value):
+        output = CCDData(self.data * value, meta=self.meta.copy(),
+                         name=self.name, uncertainty=self.uncertainty * value,
+                         mask=self.mask.copy(), memmap=self.memmap)
+        output.meta['SATURATE'] *= value
+        output.meta['GAIN'] /= value
+        output.meta['MAXLIN'] *= value
+        return output
+
     def __itruediv__(self, value):
         if isinstance(value, CCDData):
             self.uncertainty = np.abs(self.data / value.data) * \
@@ -176,7 +186,7 @@ def to_fits(self, context):
 
     def __del__(self):
         super().__del__()
-        del self.uncertainty
+        del self._uncertainty
 
     def __isub__(self, value):
         if isinstance(value, CCDData):
@@ -192,9 +202,14 @@ def __sub__(self, other):
         return type(self)(data=self.data - other.data, meta=self.meta, mask=self.mask | other.mask,
                           uncertainty=uncertainty)
 
-    def add_uncertainty(self, readnoise: np.array):
-        self._validate_array(readnoise)
-        self.uncertainty = self._init_array(readnoise)
+    @property
+    def uncertainty(self):
+        return self._uncertainty
+
+    @uncertainty.setter
+    def uncertainty(self, value: np.array):
+        self._validate_array(value)
+        self._uncertainty = self._init_array(value)
 
     def signal_to_noise(self):
         return np.abs(self.data) / self.uncertainty
@@ -302,6 +317,17 @@ def data_section(self, section):
         self.meta['DATASEC'] = section.to_region_keyword()
         self._data_section = section
 
+    @property
+    def n_sub_exposures(self):
+        n_exposures = self.meta.get('NSUBREAD', 1)
+        if str(n_exposures).lower() in ['n/a', 'unknown', 'none', '']:
+            n_exposures = 1
+        return n_exposures
+
+    @n_sub_exposures.setter
+    def n_sub_exposures(self, value):
+        self.meta['NSUBREAD'] = value
+
     def rebin(self, binning):
         # TODO: Implement me
         return self

banzai/frames.py

@@ -78,6 +78,14 @@ def background(self):
     def background(self, value):
         self.primary_hdu.background = value
 
+    @property
+    def n_sub_exposures(self):
+        return self.primary_hdu.n_sub_exposures
+
+    @n_sub_exposures.setter
+    def n_sub_exposures(self, value):
+        self.primary_hdu.n_sub_exposures = value
+
     @abc.abstractmethod
     def save_processing_metadata(self, context):
         pass
@@ -195,7 +203,10 @@ def binning(self):
         return [min(x_binnings), min(y_binnings)]
 
     def __sub__(self, other):
-        return self.primary_hdu - other.primary_hdu
+        if isinstance(other, ObservationFrame):
+            return self.primary_hdu - other.primary_hdu
+        else:
+            return self.primary_hdu - other
 
     def __isub__(self, other):
         if isinstance(other, ObservationFrame):
@@ -204,6 +215,12 @@ def __isub__(self, other):
             self.primary_hdu.__isub__(other)
         return self
 
+    def __mul__(self, other):
+        if isinstance(other, ObservationFrame):
+            return self.primary_hdu.__mul__(other.primary_hdu)
+        else:
+            return self.primary_hdu.__mul__(other)
+
     def __imul__(self, other):
         if isinstance(other, ObservationFrame):
             self.primary_hdu.__imul__(other.primary_hdu)

banzai/lco.py

@@ -450,7 +450,7 @@ def open(self, file_info, runtime_context) -> Optional[ObservationFrame]:
                     # update datasec/trimsec for fs01
                     if hdu.header.get('INSTRUME') == 'fs01':
                         self._update_fs01_sections(hdu)
-                    if hdu.data.dtype == np.uint16:
+                    if hdu.data.dtype == np.uint16 or hdu.data.dtype == np.uint32:
                         hdu.data = hdu.data.astype(np.float64)
                     # check if we need to propagate any header keywords from the primary header
                     if primary_hdu is not None:

banzai/readnoise.py

@@ -1,5 +1,6 @@
 from banzai.calibrations import CalibrationUser
 from banzai.logs import format_exception, get_logger
+import numpy as np
 
 logger = get_logger()
 
@@ -8,7 +9,7 @@ class ReadNoiseLoader(CalibrationUser):
     def apply_master_calibration(self, image, master_calibration_image):
         try:
             for image_extension, readnoise_extension in zip(image.ccd_hdus, master_calibration_image.ccd_hdus):
-                image_extension.add_uncertainty(readnoise_extension.data)
+                image_extension.uncertainty = readnoise_extension.data * np.sqrt(image_extension.n_sub_exposures)
         except:
             logger.error(f"Can't add READNOISE to image, stopping reduction: {format_exception()}", image=image)
             return None

banzai/settings.py

@@ -54,8 +54,8 @@
 EXTRA_STAGES = {'BIAS': ['banzai.bias.BiasMasterLevelSubtractor', 'banzai.bias.BiasComparer'],
                 'DARK': ['banzai.dark.DarkNormalizer', 'banzai.dark.DarkTemperatureChecker',
                          'banzai.dark.DarkComparer'],
-                'SKYFLAT': ['banzai.flats.FlatSNRChecker', 'banzai.flats.FlatNormalizer', 'banzai.qc.PatternNoiseDetector',
-                            'banzai.flats.FlatComparer'],
+                'SKYFLAT': ['banzai.flats.FlatSNRChecker', 'banzai.flats.FlatNormalizer',
+                            'banzai.qc.PatternNoiseDetector', 'banzai.flats.FlatComparer'],
                 'STANDARD': None,
                 'EXPOSE': None,
                 'EXPERIMENTAL': None}
@@ -80,7 +80,7 @@
                                   'elp': {'minute': 0, 'hour': 23},
                                   'ogg': {'minute': 0, 'hour': 3}}
 
-ASTROMETRY_SERVICE_URL = os.getenv('ASTROMETRY_SERVICE_URL', ' ')
+ASTROMETRY_SERVICE_URL = os.getenv('ASTROMETRY_SERVICE_URL', 'http://astrometry.lco.gtn/catalog')
 
 CALIBRATION_FILENAME_FUNCTIONS = {'BIAS': ('banzai.utils.file_utils.config_to_filename',
                                            'banzai.utils.file_utils.ccdsum_to_filename'),