Add scaled fwhm quantities.

Author: erykoff

decasu/configuration.py

@@ -72,6 +72,7 @@ class Configuration:
     band_field: str = 'band'
     mjd_field: str = 'mjd_obs'
     skyvar_field: str = 'skyvar'
+    fwhm_field: str = 'fwhm'
     bad_amps: Dict[int, list] = field(default_factory=_default_bad_amps)
     bad_ccds: List[int] = field(default_factory=_default_bad_ccds)
     latitude: float = -30.1690

decasu/lsst_wcs_consdb.py

@@ -7,6 +7,7 @@
 from astropy.coordinates import EarthLocation
 
 from . import decasu_globals
+from .utils import compute_visit_iqr_and_optics_scale
 
 try:
     import lsst.obs.lsst
@@ -97,8 +98,12 @@ def __init__(self, config, dbstring, bands, compute_pixels=True):
         print(f"Found {len(db_table)} detector visits for {len(self.bands)} bands.")
 
         # Add extra columns.
+        # Units of degrees.
         db_table["decasu_lst"] = np.zeros(len(db_table))
+        # Units of electrons.
         db_table["skyvar"] = db_table["sky_noise"]**2.
+        # Units of arcsec.
+        db_table[config.fwhm_field] = 2.355*config.argsec_per_pix*db_table["psf_sigma"]
 
         print("Computing local sidereal time...")
         loc = EarthLocation(lat=config.latitude*units.degree,
@@ -107,7 +112,14 @@ def __init__(self, config, dbstring, bands, compute_pixels=True):
 
         t = Time(db_table[config.mjd_field], format="mjd", location=loc)
         lst = t.sidereal_time("apparent")
-        db_table["decasu_lst"] = lst.to_value(units.degree)
+        db_table["decasu_lst"][:] = lst.to_value(units.degree)
+
+        # Compute a couple of additional psf quantities.
+        db_table[f"{config.fwhm_field}_iqr"] = np.zeros(len(db_table))
+        db_table[f"{config.fwhm_field}_optics_scale"] = np.zeros(len(db_table))
+
+        print('Computing fwhm scaled properties...')
+        compute_visit_iqr_and_optics_scale(self.config, db_table)
 
         instrument = lsst.obs.lsst.LsstCam()
         camera = instrument.getCamera()

decasu/utils.py

@@ -166,3 +166,28 @@ def read_maskfiles(expnums, maskfiles, exp_field):
             masktable = np.append(masktable, subtable[b])
 
     return masktable
+
+
+def compute_visit_iqr_and_optics_scale(config, table):
+    """Compute the visit IQR and optics scale.
+
+    Parameters
+    ----------
+    config : `decasu.Configuration`
+    table : `np.ndarray` or `astropy.table.Table`
+    """
+    u, inv = np.unique(table[config.exp_field], return_inverse=True)
+    h, rev = esutil.stat.histogram(inv, rev=True)
+    inds, = np.where(h > 0)
+
+    table[f"{config.fwhm_field}_iqr"][:] = 0.0
+    table[f"{config.fwhm_field}_optics_scale"][:] = 1.0
+
+    for ind in inds:
+        i1a = rev[rev[ind]: rev[ind + 1]]
+        isfinite = np.isfinite(table[config.fwhm_field][i1a]) & (table[config.fwhm_field][i1a] > 0.0)
+        if np.any(isfinite):
+            five, lo, hi = np.percentile(table[config.fwhm_field][i1a][isfinite], [5.0, 25.0, 75.0])
+            table[f"{config.fwhm_field}_iqr"][i1a][isfinite] = hi - lo
+
+            table[f"{config.fwhm_field}_optics_scale"][i1a] = table[f"{config.fwhm_field}"][i1a] / five

decasu/wcs_table.py

@@ -8,6 +8,7 @@
 from astropy.coordinates import EarthLocation
 
 from . import decasu_globals
+from .utils import compute_visit_iqr_and_optics_scale
 
 
 class WcsTableBuilder:
@@ -48,6 +49,8 @@ def __init__(self, config, infiles, bands, pfw_attempt_ids=[], compute_pixels=Tr
 
             # And add in the hour angle and parallactic angle fields
             dtype.extend([('decasu_lst', 'f8')])
+            dtype.extend([(f'{config.fwhm_field}_iqr', 'f4')])
+            dtype.extend([(f'{config.fwhm_field}_optics_scale', 'f4')])
 
             table = np.zeros(table_in.size, dtype=dtype)
             for name in table_in.dtype.names:
@@ -90,6 +93,9 @@ def __init__(self, config, infiles, bands, pfw_attempt_ids=[], compute_pixels=Tr
 
         print('Found %d CCDs for %d bands.' % (len(fulltable), len(bands)))
 
+        print('Computing fwhm scaled properties...')
+        compute_visit_iqr_and_optics_scale(self.config, fulltable)
+
         print('Computing local sidereal time...')
         loc = EarthLocation(lat=config.latitude*units.degree,
                             lon=config.longitude*units.degree,