Merge pull request #67 from emirkmo/refactor_photometry

Refactor photometry

flows/__init__.py

@@ -1,7 +1,7 @@
 """
 FLOWS pipeline package
 """
-from .photometry import photometry
+from .photometry import do_phot as photometry
 from .catalogs import download_catalog
 from .visibility import visibility
 from .version import get_version

flows/epsfbuilder/__init__.py

@@ -2,4 +2,4 @@
 # -*- coding: utf-8 -*-
 # flake8: noqa
 
-from .epsfbuilder import FlowsEPSFBuilder
+from .epsfbuilder import FlowsEPSFBuilder, verify_epsf

flows/epsfbuilder/epsfbuilder.py

@@ -7,8 +7,12 @@
 """
 import time
 import numpy as np
-from scipy.interpolate import griddata
+from scipy.interpolate import griddata, UnivariateSpline
 import photutils.psf
+from typing import List, Tuple
+import logging
+
+logger = logging.getLogger(__name__)
 
 
 class FlowsEPSFBuilder(photutils.psf.EPSFBuilder):
@@ -51,3 +55,31 @@ def __call__(self, *args, **kwargs):
         epsf.fit_info = dict(n_iter=len(self._epsf), max_iters=self.maxiters, time=time.time() - t0, )
 
         return epsf, stars
+
+
+def verify_epsf(epsf: photutils.psf.EPSFBuilder) -> Tuple[bool, List[float]]:
+    fwhms = []
+    epsf_ok = True
+    for a in (0, 1):
+        # Collapse the PDF along this axis:
+        profile = epsf.data.sum(axis=a)
+        itop = profile.argmax()
+        poffset = profile[itop] / 2
+
+        # Run a spline through the points, but subtract half of the peak value, and find the roots:
+        # We have to use a cubic spline, since roots() is not supported for other splines
+        # for some reason
+        profile_intp = UnivariateSpline(np.arange(0, len(profile)), profile - poffset, k=3, s=0, ext=3)
+        lr = profile_intp.roots()
+
+        # Do some sanity checks on the ePSF:
+        # It should pass 50% exactly twice and have the maximum inside that region.
+        # I.e. it should be a single gaussian-like peak
+        if len(lr) != 2 or itop < lr[0] or itop > lr[1]:
+            logger.error(f"EPSF is not a single gaussian-like peak along axis {a}")
+            epsf_ok = False
+        else:
+            axis_fwhm = lr[1] - lr[0]
+            fwhms.append(axis_fwhm)
+
+    return epsf_ok, fwhms

flows/fileio.py

@@ -0,0 +1,80 @@
+import os
+import logging
+from typing import Optional, Protocol
+from configparser import ConfigParser
+
+logger = logging.getLogger(__name__)
+
+
+class DirectoryProtocol(Protocol):
+    archive_local: str
+    output_folder: str
+
+    def set_output_dirs(self, target_name: str, fileid: int) -> None:
+        ...
+
+    def image_path(self, image_path: str) -> str:
+        ...
+
+    @property
+    def photometry_path(self) -> str:
+        ...
+
+
+class Directories:
+    """
+    Class for creating input and output directories, given a configparser instance.
+    Overwrite archive_local or output_folder to manually place elsewhere.
+    """
+    archive_local: Optional[str] = None
+    output_folder: Optional[str] = None
+
+    def __init__(self, config: ConfigParser):
+        self.config = config
+
+    def set_output_dirs(self, target_name: str, fileid: int) -> None:
+        """
+        The function is meant to be called from within a context where a
+        target_name and fileid are defined, so that the output_folder
+        can be created appropriately.
+
+        Parameters:
+            target_name (str): Target name.
+            fileid (int): The fileid of the file being processed
+        """
+
+        # Checking for None allows manual declarations to not be overwritten.
+        if self.archive_local is None:
+            self.archive_local = self._set_archive()
+
+        self.output_folder = self._set_output(target_name, fileid)
+
+        # Create output folder if necessary.
+        os.makedirs(self.output_folder, exist_ok=True)
+        logger.info("Placing output in '%s'", self.output_folder)
+
+    def _set_archive(self) -> str:
+        archive_local = self.config.get('photometry', 'archive_local', fallback=None)
+        if archive_local is not None and not os.path.isdir(archive_local):
+            raise FileNotFoundError("ARCHIVE is not available: " + archive_local)
+        logger.info(f"Using data from: {archive_local}.")
+        return archive_local
+
+    def _set_output(self, target_name: str, fileid: int) -> str:
+        """
+        Directory for output, defaults to current
+        directory if config is invalid or empty.
+        """
+        output_folder_root = self.config.get('photometry', 'output', fallback='.')
+        output_folder = os.path.join(output_folder_root, target_name, f'{fileid:05d}')
+        return output_folder
+
+    def image_path(self, image_path: str) -> str:
+        return os.path.join(self.archive_local, image_path)
+
+    @property
+    def photometry_path(self) -> str:
+        return os.path.join(self.output_folder, 'photometry.ecsv')
+
+    def save_as(self, filename: str) -> str:
+        return os.path.join(self.output_folder, filename)

flows/filters.py

@@ -0,0 +1,33 @@
+import logging
+
+FILTERS = {
+    'up': 'u_mag',
+    'gp': 'g_mag',
+    'rp': 'r_mag',
+    'ip': 'i_mag',
+    'zp': 'z_mag',
+    'B': 'B_mag',
+    'V': 'V_mag',
+    'J': 'J_mag',
+    'H': 'H_mag',
+    'K': 'K_mag',
+}
+
+FALLBACK_FILTER = 'gp'
+logger = logging.getLogger(__name__)
+
+
+def get_reference_filter(photfilter: str) -> str:
+    """
+    Translate photometric filter into table column.
+
+    Parameters:
+        photfilter (str): photometric filter corresponding to key of FILTERS
+    """
+
+    _ref_filter = FILTERS.get(photfilter, None)
+    if _ref_filter is None:
+        logger.warning(f"Could not find filter {photfilter} in catalogs. "
+                       f"Using default {FALLBACK_FILTER} filter.")
+        _ref_filter = FILTERS[FALLBACK_FILTER]
+    return _ref_filter

flows/load_image.py

@@ -10,7 +10,7 @@
 from astropy.io import fits
 from astropy.time import Time
 from astropy.wcs import WCS, FITSFixedWarning
-from typing import Tuple
+from typing import Tuple, Union, Dict, Any, Optional
 from tendrils import api
 
 from dataclasses import dataclass  # , field
@@ -20,21 +20,39 @@
 logger = logging.getLogger(__name__)  # Singleton logger instance
 
 
-# --------------------------------------------------------------------------------------------------
+@dataclass
+class InstrumentDefaults:
+    """
+    Default radius and FWHM for an instrument in arcseconds.
+    """
+    radius: float = 10
+    fwhm: float = 6.0   # Best initial guess
+    fwhm_min: float = 3.5
+    fwhm_max: float = 18.0
+
+
 @dataclass
 class FlowsImage:
     image: np.ndarray
     header: typing.Dict
-    mask: np.ndarray = None
-    peakmax: float = None
-    exptime: float = None
-    clean: np.ma.MaskedArray = None
+    mask: Optional[np.ndarray] = None
+    peakmax: Optional[float] = None
+    exptime: Optional[float] = None
+    instrument_defaults: Optional[InstrumentDefaults] = None
+    site: Optional[Dict[str, Any]] = None
+    obstime: Optional[Time] = None
+    photfilter: Optional[str] = None
+    wcs: Optional[WCS] = None
+
+    clean: Optional[np.ma.MaskedArray] = None
+    subclean: Optional[np.ma.MaskedArray] = None
+    error: Optional[np.ma.MaskedArray] = None
 
     def __post_init__(self):
         self.shape = self.image.shape
         self.wcs = self.create_wcs()
         # Make empty mask
-        if not self.mask:
+        if self.mask is None:
             self.mask = np.zeros_like(self.image, dtype='bool')
         self.check_finite()
 
@@ -169,23 +187,29 @@ class Instrument(AbstractInstrument):
     def __init__(self, image: FlowsImage = None):
         self.image = image
 
-    def get_site(self):
+    def get_site(self) -> Dict[str, Any]:
         if self.siteid is not None:
             return api.get_site(self.siteid)
 
-    def get_exptime(self):
+    def get_exptime(self) -> Union[float, int, str]:
         exptime = self.image.header.get('EXPTIME', None)
         if exptime is None:
             raise ValueError("Image exposure time could not be extracted")
         return exptime
 
-    def get_obstime(self):
+    def get_obstime(self) -> Time:
         """Default for JD, jd, utc."""
         return Time(self.image.header['JD'], format='jd', scale='utc', location=self.image.site['EarthLocation'])
 
     def get_photfilter(self):
         return self.image.header['FILTER']
 
+    def set_instrument_defaults(self):
+        """
+        Set default values for instrument.
+        """
+        self.image.instrument_defaults = InstrumentDefaults()
+
     def _get_clean_image(self):
         self.image.peakmax = self.peakmax
         self.image.site = self.get_site()
@@ -194,14 +218,16 @@ def _get_clean_image(self):
         self.image.photfilter = self.get_photfilter()
         self.image.create_masked_image()
 
-    def process_image(self, image: FlowsImage = None):
+    def process_image(self, image: FlowsImage = None) -> FlowsImage:
         """Process existing or new image."""
         if image is not None:
             self.image = image
         if self.image is None:
             raise AttributeError('No FlowsImage to be processed. Self.image was None')
 
         self._get_clean_image()
+        self.set_instrument_defaults()
+        return self.image
 
 
 class LCOGT(Instrument):
@@ -373,6 +399,14 @@ def get_photfilter(self):
                                                                         self.image.header['FILTER'])
         return photfilter
 
+class Swope_newheader(Swope):
+
+    def get_obstime(self):
+        obstime = Time(self.image.header['MJD-OBS'], format='mjd', scale='utc',
+                       location=self.image.site['EarthLocation'])
+        obstime += 0.5 * self.image.exptime * u.second
+        return obstime
+
 
 class Dupont(Instrument):
     siteid = 14
@@ -522,10 +556,35 @@ def get_obstime(self):
         return obstime
 
 
+class RATIR(Instrument):
+    siteid = 23
+
+    def get_obstime(self):
+        obstime =  Time(self.image.header['DATE-OBS'], format='isot', scale='utc',
+                        location=self.image.site['EarthLocation'])
+        return obstime
+
+
 instruments = {'LCOGT': LCOGT, 'HAWKI': HAWKI, 'ALFOSC': ALFOSC, 'NOTCAM': NOTCAM, 'PS1': PS1, 'Liverpool': Liverpool,
-               'Omega2000': Omega2000, 'Swope': Swope, 'Dupont': Dupont, 'Retrocam': RetroCam, 'Baade': Baade,
+               'Omega2000': Omega2000, 'Swope': Swope, 'Swope_newheader':Swope_newheader, 'Dupont': Dupont, 'Retrocam':
+                   RetroCam, 'Baade': Baade,
                'Sofi': Sofi, 'EFOSC': EFOSC, 'AstroNIRCam': AstroNIRCam, 'OmegaCam': OmegaCam, 'AndiCam': AndiCam,
-               'PairTel': PairTel, 'TJO': TJO}
+               'PairTel': PairTel, 'TJO': TJO, 'RATIR': RATIR, }
+
+
+def correct_barycentric(obstime: Time, target_coord: coords.SkyCoord) -> Time:
+    """
+    BARYCENTRIC CORRECTION OF TIME
+
+    Parameters:
+        obstime (astropy.time.Time): Midpoint observed image time.
+        target_coord (astropy.coords.SkyCoord): Coordinates of target in image.
+
+    Returns:
+        obstime (astropy.time.Time): Time corrected to barycenter with jpl ephemeris.
+    """
+    ltt_bary = obstime.light_travel_time(target_coord, ephemeris='jpl')
+    return obstime.tdb + ltt_bary
 
 
 def load_image(filename: str, target_coord: typing.Union[coords.SkyCoord, typing.Tuple[float, float]] = None):
@@ -602,10 +661,14 @@ def load_image(filename: str, target_coord: typing.Union[coords.SkyCoord, typing
                 instrument_name = Omega2000()
                 break
 
-            elif telescope == 'SWO' and hdr.get('SITENAME') == 'LCO':
+            elif telescope.upper().startswith('SWO') and hdr.get('SITENAME') == 'LCO':
                 instrument_name = Swope()
                 break
 
+            elif telescope.upper().startswith('SWO') and origin == 'ziggy':
+                instrument_name = Swope_newheader()
+                break
+
             elif telescope == 'DUP' and hdr.get('SITENAME') == 'LCO' and instrument == 'Direct/SITe2K-1':
                 instrument_name = Dupont()
                 break
@@ -649,10 +712,16 @@ def load_image(filename: str, target_coord: typing.Union[coords.SkyCoord, typing
                 instrument_name = TJO()
                 break
 
+            elif telescope == "OAN/SPM Harold L. Johnson 1.5-meter":
+                instrument_name = RATIR()
+                break
+
             else:
                 raise RuntimeError("Could not determine origin of image")
 
         image = FlowsImage(image=np.asarray(hdul[ext].data, dtype='float64'), header=hdr, mask=mask)
-        instrument_name.process_image(image)
+        clean_image = instrument_name.process_image(image)
 
-        return instrument_name.image
+        if target_coord is not None:
+            clean_image.obstime = correct_barycentric(clean_image.obstime, target_coord)
+        return clean_image