v0.10.8 release: Fix identified instrument bugs and rescaling improvements

VERSION

@@ -1 +1 @@
-master-v0.10.7
+master-v0.10.8

flows/instruments/__init__.py

@@ -1 +1 @@
-from .instruments import INSTRUMENTS, Instrument
+from .instruments import INSTRUMENTS, Instrument, verify_coordinates

flows/instruments/instruments.py

@@ -1,5 +1,32 @@
 """
-Load image code.
+Instrument classes that inherit from the base class.
+Modify to add your own instrument.
+
+Identifying relevant image properties:
+
+`site` = required, add manually.
+`peakmax` = optional, add manually or provide header.
+The rest of `exptime`, `obstime`, `photfilter` use the
+following (overrideable) base functions. Override
+if the default one from the baseclass does not fit
+your instrument. See:
+    ```
+    self.image.peakmax = self.peakmax
+    self.image.site = self.get_site() 
+    self.image.exptime = self.get_exptime()
+    self.image.obstime = self.get_obstime()
+    self.image.photfilter = self.get_photfilter()
+    ```
+
+Identifying the instrument for an image:
+
+Each instrument can define (one or many) of `origin`, 
+`telescope`, `instrument` fields correspinding to the 
+standard fits headers to help uniquely identify itself.
+More advanced logic is possible using `unique_headers`
+field as a dict of key,value pairs in the header. Ex:
+unique_headers = {'PRODCATG': 'SCIENCE.MEFIMAGE'}.
+These are all optional, defaults are set in baseclass.
 """
 # Standard lib
 from __future__ import annotations
@@ -72,7 +99,7 @@ class HAWKI(Instrument):
     telescope = 'ESO-VLT-U4'  # Fits Header name of TELESCOP
     instrument = 'HAWKI'  # Fits Header name of Instrument (can be partial)
     origin = 'ESO-PARANAL'  # Fits Header value of ORIGIN (if relevant)
-    unique_headers = {'PRODCATG': 'SCIENCE.MEFIMAGE'}
+    #unique_headers = {'PRODCATG': 'SCIENCE.MEFIMAGE'}
 
     def __init__(self, image: FlowsImage = None):
         super().__init__(image)
@@ -101,6 +128,10 @@ def get_ext(hdul: fits.HDUList, target_coords: coords.SkyCoord = None,
         if target_coord is None:
             raise ValueError("TARGET_COORD is needed for HAWKI images to find the correct extension")
 
+        # Incase this is not a multi-extension imageL
+        if len(hdul) == 1:
+            return 0
+
         # For HAWKI multi-extension images we search the extensions for which one contains
         # the target, Create Image from that extension.
         target_radec = [[target_coord.icrs.ra.deg, target_coord.icrs.dec.deg]]
@@ -367,7 +398,7 @@ def get_obstime(self):
 
     def get_photfilter(self):
         hdr = self.image.header
-        photfilter = {'H_Open': 'H', 'K_Open': 'K', }.get(hdr['FILTER'], hdr['FILTER'])
+        photfilter = {'H_Open': 'H', 'K_Open': 'K', 'J_Open': 'J'}.get(hdr['FILTER'], hdr['FILTER'])
         return photfilter
 
 
@@ -488,12 +519,9 @@ def get_photfilter(self):
 class Schmidt(Instrument):
     siteid = 26
     peakmax = 56_000
-    telescope = '67/91 Schmidt Telescope'  # Fits Header name of TELESCOP
-    instrument = 'Moravian G4-16000LC'  # Fits Header name of Instrument (can be partial)
-    origin = ''  # Fits Header value of ORIGIN (if relevant)
-    unique_headers = {
-        'SITELAT': 45.8494444
-    }  # Unique key value pairs from header for identifying instrument.
+    telescope = '67/91 Schmidt Telescope' 
+    instrument = 'Moravian G4-16000LC'
+    origin = '' 
 
     def get_obstime(self):
         obstime = Time(self.image.header['DATE-OBS'], format='isot', scale='utc',

flows/load_image.py

@@ -8,7 +8,7 @@
 from astropy.io import fits
 from astropy.time import Time
 
-from .instruments import INSTRUMENTS
+from .instruments import INSTRUMENTS, verify_coordinates
 from .image import FlowsImage
 from .utilities import create_logger
 logger = create_logger()
@@ -39,6 +39,7 @@ def load_image(filename: str, target_coord: Union[coords.SkyCoord, Tuple[float,
         for inst_name, inst_cls in INSTRUMENTS:
             if inst_cls.identifier(telescope, origin, instrument, hdr):
                 logger.info(f"Image is using instrument {inst_name}")
+                target_coord = verify_coordinates(target_coord)
                 ext = inst_cls.get_ext(hdul, target_coord)
                 mask = inst_cls.get_mask(hdul)
                 # Default = None is to only mask all non-finite values, override here is additive.

flows/photometry.py

@@ -316,10 +316,11 @@ def diff_psf_phot(self) -> Table:
                                              init_table=self.init_guesses_diff.init_guess_diff)
         return psf_tbl
 
-    def rescale_uncertainty(self, psfphot_tbl: Table, dynamic: bool = True, static_fwhm: float = 2.5):
+    def rescale_uncertainty(self, psfphot_tbl: Table, dynamic: bool = True, 
+                            static_fwhm: float = 2.5, epsilon_mag: float = 0.004,
+                            ensure_greater: bool = True):
         """
-        Rescale the uncertainty of the PSF photometry to match the uncertainty of the
-        photometry.
+        Rescale the uncertainty of the PSF photometry using a variable fitsize.
 
         Parameters
         ----------
@@ -329,6 +330,9 @@ def rescale_uncertainty(self, psfphot_tbl: Table, dynamic: bool = True, static_f
             Dynamically decide FWHM multiple for rescaling.
         static_fwhm : float
             FWHM multiple to use incase dynamic fails or don't want to use it. Default 2.5 determined empirically.
+        epsilon_mag : float
+            Small magnitude change within which new and old uncertainties are considered the same. 
+            Should be smaller than ~1/2 the expected uncertainty.
         """
         # Rescale psf errors from fit iteratively
         fit_shapes = self.photometry.get_fit_shapes(self.fwhm, self.psf_builder.star_size)
@@ -361,6 +365,12 @@ def rescale_uncertainty(self, psfphot_tbl: Table, dynamic: bool = True, static_f
         logger.info(f"Recalculating all reference uncertainties using new fitsize:"
                     f" {fit_shape} pixels, ({fit_shape/self.fwhm if dynamic else static_fwhm :.2} * FWHM).")
         psfphot_tbl_rescaled = self.psfphot(fit_shape)
+        if psfphot_tbl['flux_unc'][0] > psfphot_tbl_rescaled['flux_unc'][0] + epsilon_mag and ensure_greater:
+            logger.info("Recalculated uncertainties were smaller than original and ``ensure_greater`` was True:"
+                        "Not using rescaled uncertainties for the SN.")
+            psfphot_tbl['flux_unc'][1:] = psfphot_tbl_rescaled['flux_unc'][1:]
+            return psfphot_tbl
+
         psfphot_tbl['flux_unc'] = psfphot_tbl_rescaled['flux_unc']
         return psfphot_tbl
 
@@ -415,7 +425,7 @@ def create_from_fid(cls, fid: int, directories: Optional[DirectoryProtocol] = No
 
 def do_phot(fileid: int, cm_timeout: Optional[float] = None, make_plots: bool = True,
             directories: Optional[DirectoryProtocol] = None, datafile: Optional[Dict[str, Any]] = None,
-            rescale_dynamic: bool = True) -> ResultsTable:
+            rescale: bool = True, rescale_dynamic: bool = True) -> ResultsTable:
     # Set up photometry runner
     pm = PhotometryManager.create_from_fid(fileid, directories=directories, datafile=datafile, create_directories=True)
 
@@ -430,8 +440,10 @@ def do_phot(fileid: int, cm_timeout: Optional[float] = None, make_plots: bool =
 
     # Do photometry
     apphot_tbl = pm.apphot()
-    psfphot_tbl = ResultsTable.verify_uncertainty_column(pm.psfphot())  # Verify uncertainty exists after PSF phot.
-    psfphot_tbl = pm.rescale_uncertainty(psfphot_tbl, dynamic=rescale_dynamic)  # Rescale uncertainties
+    # Verify uncertainty exists after PSF phot:
+    psfphot_tbl = ResultsTable.verify_uncertainty_column(pm.psfphot())  
+    if rescale:  # Rescale uncertainties
+        psfphot_tbl = pm.rescale_uncertainty(psfphot_tbl, dynamic=rescale_dynamic) 
 
     # Build results table and calculate magnitudes
     pm.make_result_table(psfphot_tbl, apphot_tbl)
@@ -445,10 +457,11 @@ def do_phot(fileid: int, cm_timeout: Optional[float] = None, make_plots: bool =
 
 def timed_photometry(fileid: int, cm_timeout: Optional[float] = None, make_plots: bool = True,
                      directories: Optional[DirectoryProtocol] = None, save: bool = True,
-                     datafile: Optional[Dict[str, Any]] = None, rescale_dynamic: bool = True) -> ResultsTable:
+                     datafile: Optional[Dict[str, Any]] = None, rescale: bool = True,
+                     rescale_dynamic: bool = True) -> ResultsTable:
     # TODO: Timer should be moved out of this function.
     tic = default_timer()
-    results_table = do_phot(fileid, cm_timeout, make_plots, directories, datafile, rescale_dynamic)
+    results_table = do_phot(fileid, cm_timeout, make_plots, directories, datafile, rescale, rescale_dynamic)
 
     # Save the results table:
     if save:

run_photometry.py

@@ -11,8 +11,8 @@
 from flows.utilities import create_logger, parse_log_level, create_warning_logger, remove_file_handlers
 
 
-def process_fileid(fid, autoupload=False, cm_timeout=None, no_plots=False,
-                   rescale_dynamic=True) -> result_model.ResultsTable:
+def process_fileid(fid, autoupload: bool = False, cm_timeout=None, no_plots: bool = False,
+                   rescale: bool = True, rescale_dynamic: bool = True) -> result_model.ResultsTable:
     # Create the output directory if it doesn't exist:
     datafile = api.get_datafile(fid)
     directories = fileio.Directories.from_fid(fid, datafile=datafile)
@@ -25,7 +25,8 @@ def process_fileid(fid, autoupload=False, cm_timeout=None, no_plots=False,
             api.set_photometry_status(fid, 'running')
 
         table = photometry(fileid=fid, cm_timeout=cm_timeout, make_plots=not no_plots,
-                           directories=directories, datafile=datafile, rescale_dynamic=rescale_dynamic)
+                           directories=directories, datafile=datafile, rescale=rescale, 
+                           rescale_dynamic=rescale_dynamic)
 
     except (SystemExit, KeyboardInterrupt):
         logger.error("Aborted by user or system.")
@@ -79,6 +80,8 @@ def main():
     group.add_argument('--fixposdiff', action='store_true',
                        help="Fix SN position during PSF photometry of difference image. "
                             "Useful when difference image is noisy.")
+    group.add_argument('--rescale-off', action='store_false',
+                        help='Turn off uncertainty rescaling.')
     group.add_argument('--wcstimeout', type=int, default=None, help="Timeout in Seconds for WCS.")
     args = parser.parse_args()
 
@@ -118,8 +121,9 @@ def main():
     fileids = list(set(fileids))
 
     # Create function wrapper:
-    process_fileid_wrapper = functools.partial(process_fileid, autoupload=args.autoupload, cm_timeout=args.wcstimeout,
-                                               no_plots=args.noplots, rescale_dynamic=not args.rescale_static)
+    process_fileid_wrapper = functools.partial(process_fileid, autoupload=args.autoupload,
+                                               cm_timeout=args.wcstimeout, no_plots=args.noplots,
+                                               rescale=args.rescale_off, rescale_dynamic=not args.rescale_static)
 
     if threads > 1:
         # process in parallel: