Add weight and flags image creation to preprocessing

default_config.yaml

@@ -111,3 +111,4 @@ DECam:
     linearity: DECamMasterCal_56475/linearity/linearity_table_v0.4.fits
     fringebase: DECamMasterCal_56876/fringecor/DECam_Master_20131115v1
     flatbase: DECamMasterCal_56876/starflat/DECam_Master_20130829v3
+    bpmbase: DECamMasterCal_56876/bpm/DECam_Master_20140209v2_cd_

models/decam.py

@@ -24,6 +24,7 @@
 from util.config import Config
 import util.util
 from util.retrydownload import retry_download
+from models.enums_and_bitflags import string_to_bitflag, flag_image_bits_inverse
 
 class DECam(Instrument):
 
@@ -222,6 +223,35 @@ def get_ra_dec_for_section( self, exposure, section_id ):
         return ( exposure.ra + self._chip_radec_off[section_id]['dra'] / math.cos( exposure.dec * math.pi / 180. ),
                  exposure.dec + self._chip_radec_off[section_id]['ddec'] )
 
+    def get_standard_flags_image( self, section_id ):
+        # NOTE : there's a race condition here; multiple
+        # processes might try to locally cache the
+        # file all at the same time.
+        #
+        # For now, just not going to worry about it.
+        #
+        cfg = Config.get()
+        ccdnum = f'{self._chip_radec_off[section_id]["ccdnum"]:02d}'
+        rempath = pathlib.Path( f'{cfg.value("DECam.calibfiles.bpmbase")}{ccdnum}.fits' )
+        filepath = pathlib.Path( FileOnDiskMixin.local_path ) / "DECam_default_calibrators" / "bpm" / rempath.name
+
+        if not filepath.exists():
+            url = f'{cfg.value("DECam.calibfiles.urlbase")}{str(rempath)}'
+            retry_download( url, filepath )
+
+        with fits.open( filepath, memmap=False ) as hdu:
+            rawbpm = hdu[0].data
+
+        # TODO : figure out what the bits mean in this bad pixel mask file!
+        # For now, call anything non-zero as "bad"
+        # (If we ever change this, we have to fix the
+        # flag_image_bits dictionary in enums_and_bitflags,
+        # and everything that has used it...)
+        bpm = np.zeros( rawbpm.shape, dtype=np.uint16 )
+        bpm[ rawbpm != 0 ] = string_to_bitflag( 'bad pixel', flag_image_bits_inverse )
+
+        return bpm
+
     def get_gain_at_pixel( self, image, x, y, section_id=None ):
         if image is None:
             return Instrument.get_gain_at_pixel( image, x, y, section_id=section_id )
@@ -239,6 +269,19 @@ def get_gain_at_pixel( self, image, x, y, section_id=None ):
         else:
             return float( image.header[ 'GAINB' ] )
 
+    def average_gain( self, image, section_id=None ):
+        if image is None:
+            return Instrument.average_again( self, None, section_id=section_id )
+        return ( float( image.raw_header['GAINA'] ) + float( image.raw_header['GAINB'] ) ) / 2.
+
+
+    def average_saturation_limit( self, image, section_id=None ):
+        if image is None:
+            return Instrument.average_saturation_limit( self, image, section_id=section_id )
+        # Although the method name is "average...", return the lower saturation
+        #  limit to be conservative
+        return min( float( image.raw_header['SATURATA'] ), float( image.raw_header['SATURATB'] ) )
+
     @classmethod
     def _get_fits_hdu_index_from_section_id(cls, section_id):
         """

models/enums_and_bitflags.py

@@ -335,3 +335,11 @@ def string_to_bitflag(value, dictionary):
     6: 'illumination'
 }
 image_preprocessing_inverse = {EnumConverter.c(v):k for k, v in image_preprocessing_dict.items()}
+
+# bitflag used in flag images
+flag_image_bits = {
+    0: 'bad pixel',        # Bad pixel flagged by the instrument
+    1: 'zero weight',
+    2: 'saturated',
+}
+flag_image_bits_inverse = { EnumConverter.c(v):k for k, v in flag_image_bits.items() }

models/image.py

@@ -707,6 +707,14 @@ def from_exposure(cls, exposure, section_id):
         # read the header from the exposure file's individual section data
         new._raw_header = exposure.section_headers[section_id]
 
+        # Because we will later be writing out float data (BITPIX=-32)
+        # -- or whatever the type of raw_data is -- we have to make sure
+        # there aren't any vestigal BSCALE and BZERO keywords in the
+        # header.
+        for delkw in [ 'BSCALE', 'BZERO' ]:
+            if delkw in new.raw_header:
+                del new.raw_header[delkw]
+
         # numpy array axis ordering is backwards from FITS ordering
         width = new.raw_data.shape[1]
         height = new.raw_data.shape[0]

models/instrument.py

@@ -370,6 +370,8 @@ def __init__(self, **kwargs):
         self.gain = getattr(self, 'gain', None)  # gain in electrons/ADU (e.g., 4.0)
         self.saturation_limit = getattr(self, 'saturation_limit', None)  # saturation limit in electrons (e.g., 100000)
         self.non_linearity_limit = getattr(self, 'non_linearity_limit', None)  # non-linearity limit in electrons
+        self.background_box_size = getattr( self, 'background_box_size', 256 ) # Box size for sep background estimation
+        self.background_filt_size = getattr( self, 'background_filt_size', 3 ) # Filter size for sep background
 
         self.allowed_filters = getattr(self, 'allowed_filters', None)  # list of allowed filter (e.g., ['g', 'r', 'i'])
 
@@ -411,27 +413,37 @@ def __repr__(self):
 
     @classmethod
     def get_section_ids(cls):
-        """
-        Get a list of SensorSection identifiers for this instrument.
+        """Get a list of SensorSection identifiers for this instrument.
+
+        Returns
+        -------
+           list of str
 
         THIS METHOD MUST BE OVERRIDEN BY THE SUBCLASS.
         """
         raise NotImplementedError("This method must be implemented by the subclass.")
 
     @classmethod
     def check_section_id(cls, section_id):
-        """
-        Check that the type and value of the section is compatible with the instrument.
-        For example, many instruments will key the section by a running integer (e.g., CCD ID),
-        while others will use a string (e.g., channel 'A').
+        """Check that the type and value of the section is compatible with the instrument.
+
+        For example, many instruments will key the section by a running
+        integer (e.g., CCD ID), while others will use a string (e.g.,
+        channel 'A').
 
         Will raise a meaningful error if not compatible.
 
         Subclasses should override this method to be more specific
-        (e.g., to test if an integer is in range).
+        (e.g., to test if an integer is in range).  IMPORTANT NOTE:
+        subclasses must *always* be able to handle a string section_id,
+        even if that subclass only uses integer section_ids.  (Reason:
+        the "section_id" field of the Image model is type string, so
+        sometimes strings will come in.)  See DemoInstrument as an example.
 
         THIS METHOD CAN BE OVERRIDEN TO MAKE THE CHECK MORE SPECIFIC
-        (e.g., to check only for integers, to check the id is in range).
+        (e.g., to verify that the actual section_id is valid for the
+        instrument)
+
         """
         if not isinstance(section_id, (int, str)):
             raise ValueError(f"The section_id must be an integer or string. Got {type(section_id)}. ")
@@ -488,7 +500,7 @@ def get_section(self, section_id):
         if self.sections is None:
             raise RuntimeError("No sections loaded for this instrument. Use fetch_sections() first.")
 
-        return self.sections.get(section_id)
+        return self.sections.get( str(section_id) )
 
     def fetch_sections(self, session=None, dateobs=None):
         """
@@ -552,7 +564,7 @@ def fetch_sections(self, session=None, dateobs=None):
                     ).all()
 
             for sid in self.get_section_ids():
-                sec = [s for s in all_sec if s.identifier == str(sid)]
+                sec = [s for s in all_sec if s.identifier == sid]
                 if len(sec) > 0:
                     self.sections[sid] = sec[0]
                 else:
@@ -613,6 +625,7 @@ def get_property(self, section_id, prop):
         """
 
         section = self.get_section(section_id)
+        section_id = section.identifier
         if section is not None:
             if hasattr(section, prop) and getattr(section, prop) is not None:
                 return getattr(section, prop)
@@ -932,7 +945,7 @@ def get_standard_flags_image( self, section_id ):
         """
 
         sec = self.get_section( section_id )
-        return np.zeros( [ sec.size_y, sec.size_x ], dtype=numpy.uint16 )
+        return np.zeros( [ sec.size_y, sec.size_x ], dtype=np.uint16 )
 
     def get_gain_at_pixel( self, image, x, y, section_id=None ):
         """Get the gain of an image at a given pixel position.
@@ -946,23 +959,72 @@ def get_gain_at_pixel( self, image, x, y, section_id=None ):
         Parameters
         ----------
         image: Image or None
-          The Image.  If None, section_id must not be none.
+          The Image.  If None, will look at section_id
         x, y: int or float
           Position on the image in C-coordinates (0 offset).  Remember
           that numpy arrays are indexed [y, x].
         section_id:
-          If Image is None, pass a non-null section_id to get the default gain.
+          Ignored if image is not None.  If image is None, will get the
+          default gain for this section.  If both this and image are
+          None, will return the instrument default gain, or 1 if there
+          is no instrument default gain.
 
         Returns
         -------
         float
 
         """
-        sec = self.get_section( section_id if image is None else image.section_id )
-        if sec.gain is None:
+        if image is not None:
+            return self.get_section( image.section_id ).gain
+        elif section_id is not None:
+            return self.get_section( section_id ).gain
+        elif sec.gain is not None:
             return self.gain
         else:
-            return sec.gain
+            return 1.
+
+    def average_gain( self, image, section_id=None ):
+        """Get an average gain for the image.
+
+        THIS SHOULD USUALLY BE OVERRIDDEN BY SUBCLASSES.  By default,
+        it's going to assume that the gain property of the sensor
+        section is good, or if it's null, that the gain property of the
+        instrument is good.  Subclasses should use the image header
+        information.
+
+        Parameters
+        ----------
+        image: Image or None
+          The Image.  If None, section_id must not be none
+        section_id:
+          If Image is None, pass a non-null section_id to get the default gain.
+
+        Returns
+        -------
+        float
+
+        """
+        return self.get_gain_at_pixel( image, 0, 0, section_id=section_id )
+
+    def average_saturation_limit( self, image, section_id=None ):
+        """Get an average saturation limit in ADU for the image.
+
+        THIS SHOULD USUALLY BE OVERRIDDEN BY SUBCLASSES, for the same
+        reason as average_gain.
+
+        Parameters
+        ----------
+        image: Image or None
+          The Image.  If None, section_id must not be none
+        section_id: int or str
+          If Image is None, pass a non-null section_id to get the default gain.
+
+        Returns
+        -------
+        float
+
+        """
+        return self.saturation_limit
 
     @classmethod
     def _get_header_keyword_translations(cls):
@@ -1010,18 +1072,21 @@ def _get_header_values_converters(cls):
 
     @classmethod
     def _get_fits_hdu_index_from_section_id(cls, section_id):
-        """
-        Translate the section_id into the index of the HDU in the FITS file.
-        For example, if we have an instrument with 10 CCDs, numbered 0 to 9,
-        the HDU list will probably contain a generic HDU at index 0,
-        and the individual section information in 1 through 10, so
-        the function should return section_id+1.
-        Another example could have section_id=A give an index 1,
-        and a section_id=B give an index 2 (so the function will read
-        from a dictionary to translate the values).
-
-        THIS METHOD SHOULD BE OVERRIDEN BY SUBCLASSES,
-        in particular when section_id is a string.
+        """Translate the section_id into the index of the HDU in the FITS file.
+
+        For example, if we have an instrument with 10 CCDs, numbered 0
+        to 9, the HDU list will probably contain a generic HDU at index
+        0, and the individual section information in 1 through 10, so
+        the function should return section_id+1.  Another example could
+        have section_id=A give an index 1, and a section_id=B give an
+        index 2 (so the function will read from a dictionary to
+        translate the values).
+
+        THIS METHOD SHOULD USUALLY BE OVERRIDEN BY SUBCLASSES.  Only in
+        the special case where the section IDs are [0, 1, 2, 3, ...]
+        and in an exposure those correspond to FITS HDUs [1, 2, 3, 4,
+        ...] (or indexes into an astropy.io.fits HDU that is the same as
+        the section ID).  will this implementation work.
 
         Parameters
         ----------
@@ -1031,10 +1096,10 @@ def _get_fits_hdu_index_from_section_id(cls, section_id):
         Returns
         -------
         index: int
-            The index of the HDU in the FITS file.
-            Note that the index is 0-based, as this value is
-            used in the astropy.io.fits functions/objects,
-            not in the native FITS format (which is 1-based).
+            The index of the HDU in the FITS file.  Note that the index
+            is 1-based, as it corresponds to the index in the native
+            FITS file.
+
         """
         cls.check_section_id(section_id)
         return int(section_id) + 1
@@ -1588,21 +1653,20 @@ def __init__(self, **kwargs):
 
     @classmethod
     def get_section_ids(cls):
+        """Get a list of SensorSection identifiers for this instrument.
 
+        See Instrument.get_section_ids for interface.
         """
-        Get a list of SensorSection identifiers for this instrument.
-        """
-        return [0]
+
+        return [ '0' ]
 
     @classmethod
     def check_section_id(cls, section_id):
         """
         Check if the section_id is valid for this instrument.
         The demo instrument only has one section, so the section_id must be 0.
         """
-        if not isinstance(section_id, int):
-            raise ValueError(f"section_id must be an integer. Got {type(section_id)} instead.")
-        if section_id != 0:
+        if ( not isinstance(section_id, (str, int)) ) or ( int(section_id) != 0 ):
             raise ValueError(f"section_id must be 0 for this instrument. Got {section_id} instead.")
 
     def _make_new_section(self, identifier):