Wirc full pipeline (#412)

Co-authored-by: Robert Stein <rdstein@caltech.edu>

mirar/catalog/gaia.py

@@ -6,6 +6,7 @@
 
 import astropy.table
 import astropy.units as u
+import numpy as np
 from astropy.coordinates import SkyCoord
 from astroquery.gaia import Gaia
 
@@ -26,18 +27,40 @@ def __init__(
         self,
         *args,
         filter_name: str = "j",
-        ph_qual_cut: bool = False,
         snr_threshold: float = 5,
         trim: bool = False,
         image_catalog_path: Optional[str] = None,
+        acceptable_j_ph_quals: str | list[str] = None,
+        acceptable_h_ph_quals: str | list[str] = None,
+        acceptable_k_ph_quals: str | list[str] = None,
         **kwargs,
     ):
         super().__init__(*args, filter_name=filter_name, **kwargs)
-        self.ph_qual_cut = ph_qual_cut
+
         self.trim = trim
         self.image_catalog_path = image_catalog_path
         self.snr_threshold = snr_threshold
 
+        if isinstance(acceptable_j_ph_quals, str):
+            acceptable_j_ph_quals = [acceptable_j_ph_quals]
+        if isinstance(acceptable_h_ph_quals, str):
+            acceptable_h_ph_quals = [acceptable_h_ph_quals]
+        if isinstance(acceptable_k_ph_quals, str):
+            acceptable_k_ph_quals = [acceptable_k_ph_quals]
+
+        self.acceptable_ph_quals = {
+            "j": acceptable_j_ph_quals,
+            "h": acceptable_h_ph_quals,
+            "k": acceptable_k_ph_quals,
+        }
+
+        if self.acceptable_ph_quals[self.filter_name.lower()] is None:
+            self.acceptable_ph_quals[self.filter_name.lower()] = ["A"]
+
+        for filt in self.acceptable_ph_quals:
+            if self.acceptable_ph_quals[filt] is None:
+                self.acceptable_ph_quals[filt] = ["A", "B", "C"]
+
         logger.debug(f"Sextractor catalog path is {self.image_catalog_path}")
 
     def get_catalog(
@@ -62,14 +85,9 @@ def get_catalog(
             f"AND tmass.{self.filter_name}_m > {self.min_mag:.2f} "
             f"AND tmass.{self.filter_name}_m < {self.max_mag:.2f} "
             f"AND tbest.number_of_mates=0 "
-            f"AND tbest.number_of_neighbours=1"
+            f"AND tbest.number_of_neighbours=1;"
         )
 
-        if self.ph_qual_cut:
-            cmd += "AND tmass.ph_qual='AAA';"
-        else:
-            cmd += ";"
-
         job = Gaia.launch_job_async(cmd, dump_to_file=False)
         src_list = job.get_results()
         src_list["ph_qual"] = src_list["ph_qual"].astype(str)
@@ -80,6 +98,17 @@ def get_catalog(
         src_list["magnitude_err"] = src_list[f"{self.filter_name.lower()}_msigcom"]
         logger.info(f"Found {len(src_list)} sources in Gaia")
 
+        j_phquals = [x[0] for x in src_list["ph_qual"]]
+        h_phquals = [x[0] for x in src_list["ph_qual"]]
+        k_phquals = [x[0] for x in src_list["ph_qual"]]
+
+        j_phmask = np.array([x in self.acceptable_ph_quals["j"] for x in j_phquals])
+        h_phmask = np.array([x in self.acceptable_ph_quals["h"] for x in h_phquals])
+        k_phmask = np.array([x in self.acceptable_ph_quals["k"] for x in k_phquals])
+
+        phmask = j_phmask & h_phmask & k_phmask
+
+        src_list = src_list[phmask]
         src_list = src_list[src_list["magnitude_err"] < 1.086 / self.snr_threshold]
         if self.trim:
             if self.image_catalog_path is None:

mirar/io.py

@@ -11,6 +11,8 @@
 from astropy.io import fits
 from astropy.utils.exceptions import AstropyUserWarning
 
+from mirar.paths import BASE_NAME_KEY, RAW_IMG_KEY
+
 
 def create_fits(data: np.ndarray, header: fits.Header | None) -> fits.PrimaryHDU:
     """
@@ -70,6 +72,12 @@ def open_fits(path: str | Path) -> tuple[np.ndarray, fits.Header]:
         data = hdu.data
         header = hdu.header
 
+        if BASE_NAME_KEY not in header:
+            header[BASE_NAME_KEY] = Path(path).name
+
+        if RAW_IMG_KEY not in header.keys():
+            header[RAW_IMG_KEY] = path
+
     return data, header
 
 

mirar/paths.py

@@ -238,6 +238,7 @@ def get_astrometry_keys() -> list:
 UNC_IMG_KEY = "UNCPATH"
 PROC_HISTORY_KEY = "CALSTEPS"
 PROC_FAIL_KEY = "PROCFAIL"
+ASTROMETRY_FILE_KEY = "ASTRFILE"
 LATEST_SAVE_KEY = "SAVEPATH"
 LATEST_WEIGHT_SAVE_KEY = "WGHTPATH"
 SEXTRACTOR_HEADER_KEY = "SRCCAT"
@@ -252,6 +253,7 @@ def get_astrometry_keys() -> list:
 GAIN_KEY = "GAIN"
 EXPTIME_KEY = "AEXPTIME"
 ZP_KEY = "ZP"
+SATURATE_KEY = "SATURATE"
 PSF_FLUX_KEY = "psf_flux"
 PSF_FLUXUNC_KEY = "psf_fluxunc"
 MAG_PSF_KEY = "magpsf"
@@ -261,13 +263,14 @@ def get_astrometry_keys() -> list:
 CAND_NAME_KEY = "objectId"
 CAND_RA_KEY = "ra"
 CAND_DEC_KEY = "dec"
+FITS_MASK_KEY = "MASKFITS"
 sextractor_checkimg_keys = {
     "BACKGROUND": "BKGPT",
     "BACKGROUND_RMS": "BKGRMS",
     "MINIBACKGROUND": "MINIBKG",
     "MINIBACK_RMS": "MINIBGRM",
 }
-
+STACKED_COMPONENT_IMAGES_KEY = "COMPENTS"
 APFLUX_PREFIX_KEY = "fluxap"
 APFLUXUNC_PREFIX_KEY = "fluxucap"
 APMAG_PREFIX_KEY = "magap"

mirar/pipelines/__init__.py

@@ -1,3 +1,6 @@
+"""
+Central location for all pipelines. This is where you should add new pipelines.
+"""
 import logging
 
 from mirar.errors import ProcessorError
@@ -14,19 +17,30 @@
 
 
 class PipelineConfigError(ProcessorError, KeyError):
-    pass
+    """
+    Error raised when a pipeline is not found
+    """
 
 
-def get_pipeline(instrument, selected_configurations=None, *args, **kwargs):
+def get_pipeline(instrument, selected_configurations=None, *args, **kwargs) -> Pipeline:
+    """
+    Function to get pipeline
+
+    :param instrument: Name of instrument
+    :param selected_configurations: Configurations to use
+    :param args: args
+    :param kwargs: kwargs
+    :return: pipeline
+    """
     try:
         pipeline = Pipeline.pipelines[instrument.lower()]
         logger.info(f"Found {instrument} pipeline")
-    except KeyError:
+    except KeyError as exc:
         err = (
             f"Unrecognised pipeline {instrument}. "
             f"Available pipelines are: {Pipeline.pipelines.keys()}"
         )
         logger.error(err)
-        raise PipelineConfigError(err)
+        raise PipelineConfigError(err) from exc
 
     return pipeline(selected_configurations=selected_configurations, *args, **kwargs)

mirar/pipelines/sedmv2/blocks.py

@@ -23,10 +23,10 @@
 )
 from mirar.pipelines.sedmv2.load_sedmv2_image import load_raw_sedmv2_image
 from mirar.processors import BiasCalibrator, FlatCalibrator
-from mirar.processors.anet import AstrometryNet
 from mirar.processors.astromatic import PSFex, Sextractor, Swarp
+from mirar.processors.astrometry.anet import AstrometryNet
 from mirar.processors.csvlog import CSVLog
-from mirar.processors.mask import MaskPixels
+from mirar.processors.mask import MaskPixelsFromPath
 from mirar.processors.photcal import PhotCalibrator
 from mirar.processors.photometry.aperture_photometry import (
     CandidateAperturePhotometry,
@@ -79,7 +79,7 @@
 ]  # pylint: disable=duplicate-code
 
 reduce = [
-    MaskPixels(mask_path=sedmv2_mask_path),
+    MaskPixelsFromPath(mask_path=sedmv2_mask_path),
     BiasCalibrator(),
     ImageSelector(("OBSTYPE", ["FLAT", "SCIENCE"])),
     ImageBatcher(split_key="filter"),

mirar/pipelines/summer/__init__.py

@@ -1 +1,4 @@
+"""
+Pipline for SUMMER data
+"""
 from mirar.pipelines.summer.summer_pipeline import SummerPipeline

mirar/pipelines/summer/blocks.py

@@ -23,6 +23,7 @@
 from mirar.pipelines.summer.generator import (
     summer_astrometric_catalog_generator,
     summer_photometric_catalog_generator,
+    summer_photometric_img_catalog_purifier,
     summer_reference_image_generator,
     summer_reference_image_resampler,
     summer_reference_psfex,
@@ -35,7 +36,7 @@
 from mirar.pipelines.summer.models import Exposure, Proc, Raw
 from mirar.processors import BiasCalibrator, FlatCalibrator
 from mirar.processors.astromatic import PSFex, Scamp, Sextractor, Swarp
-from mirar.processors.autoastrometry import AutoAstrometry
+from mirar.processors.astrometry.autoastrometry import AutoAstrometry
 from mirar.processors.candidates.candidate_detector import DetectCandidates
 from mirar.processors.candidates.utils import DataframeWriter, RegionsWriter
 from mirar.processors.cosmic_rays import LACosmicCleaner
@@ -44,7 +45,7 @@
     DatabaseImageExporter as PSQLDatabaseImageExporter,
 )
 from mirar.processors.database.database_modifier import ModifyImageDatabaseSeq
-from mirar.processors.mask import MaskPixels
+from mirar.processors.mask import MaskPixelsFromPath
 from mirar.processors.photcal import PhotCalibrator
 from mirar.processors.photometry.aperture_photometry import CandidateAperturePhotometry
 from mirar.processors.photometry.psf_photometry import CandidatePSFPhotometry
@@ -123,7 +124,7 @@
         duplicate_protocol="replace",
         q3c_bool=False,
     ),
-    MaskPixels(mask_path=summer_mask_path),
+    MaskPixelsFromPath(mask_path=summer_mask_path),
     DatabaseImageExporter(db_table=Raw, duplicate_protocol="replace", q3c_bool=False),
     ImageSelector(("OBSTYPE", ["BIAS", "FLAT", "SCIENCE"])),
 ]
@@ -133,7 +134,7 @@
 ]
 
 test_cr = [
-    MaskPixels(mask_path=summer_mask_path),
+    MaskPixelsFromPath(mask_path=summer_mask_path),
     BiasCalibrator(),
     ImageSelector(("OBSTYPE", ["FLAT", "SCIENCE"])),
     ImageBatcher(split_key="filter"),
@@ -179,7 +180,10 @@
         checkimage_type="BACKGROUND_RMS",
         **sextractor_photometry_config
     ),
-    PhotCalibrator(ref_catalog_generator=summer_photometric_catalog_generator),
+    PhotCalibrator(
+        ref_catalog_generator=summer_photometric_catalog_generator,
+        image_photometric_catalog_purifier=summer_photometric_img_catalog_purifier,
+    ),
     ImageSaver(
         output_dir_name="processed",
         # TODO: work out why this was ever here...

mirar/pipelines/summer/config/files/astrom.sex

@@ -58,7 +58,6 @@ PHOT_PETROPARAMS 1.0,2.0                        # MAG_PETRO parameters: <Petrosi
                                 # <min_radius>
 
 SATUR_KEY        SATURATE       # keyword for saturation level (in ADUs)
-SATUR_LEVEL      60000        # level (in ADUs) at which arises saturation
 
 MAG_ZEROPOINT    0            # magnitude zero-point
 MAG_GAMMA        4.0            # gamma of emulsion (for photographic scans)

mirar/pipelines/summer/config/files/photomCat.sex

@@ -54,8 +54,8 @@ PHOT_AUTOPARAMS  2.5, 3.5       # MAG_AUTO parameters: <Kron_fact>,<min_radius>
 PHOT_PETROPARAMS 2.0, 3.5       # MAG_PETRO parameters: <Petrosian_fact>,
                                 # <min_radius>
 
-SATUR_KEY        SATURATE       # keyword for saturation level (in ADUs) 
-SATUR_LEVEL      25000        # level (in ADUs) at which arises saturation
+SATUR_KEY        SATURATE       # keyword for saturation level (in ADUs)
+SATUR_LEVEL      1000000000.0       # level (in ADUs) at which arises saturation
 
 MAG_ZEROPOINT    0            # magnitude zero-point
 MAG_GAMMA        4.0            # gamma of emulsion (for photographic scans)

mirar/pipelines/summer/generator.py

@@ -4,6 +4,8 @@
 """
 import logging
 
+from astropy.table import Table
+
 from mirar.catalog import BaseCatalog, Gaia2Mass
 from mirar.catalog.vizier import PS1, SkyMapper
 from mirar.catalog.vizier.sdss import SDSS, NotInSDSSError, in_sdss
@@ -36,10 +38,34 @@ def summer_astrometric_catalog_generator(image: Image) -> Gaia2Mass:
         trim=True,
         image_catalog_path=temp_cat_path,
         filter_name="j",
+        acceptable_j_ph_quals=["A", "B", "C"],
     )
     return cat
 
 
+def summer_photometric_img_catalog_purifier(catalog: Table, image: Image) -> Table:
+    """
+    Default function to purify the photometric image catalog
+    """
+    edge_width_pixels = 100
+    fwhm_threshold_arcsec = 4.0
+    x_lower_limit = edge_width_pixels
+    x_upper_limit = image.get_data().shape[1] - edge_width_pixels
+    y_lower_limit = edge_width_pixels
+    y_upper_limit = image.get_data().shape[0] - edge_width_pixels
+
+    clean_mask = (
+        (catalog["FLAGS"] == 0)
+        & (catalog["FWHM_WORLD"] < fwhm_threshold_arcsec / 3600.0)
+        & (catalog["X_IMAGE"] > x_lower_limit)
+        & (catalog["X_IMAGE"] < x_upper_limit)
+        & (catalog["Y_IMAGE"] > y_lower_limit)
+        & (catalog["Y_IMAGE"] < y_upper_limit)
+    )
+
+    return catalog[clean_mask]
+
+
 def summer_photometric_catalog_generator(image: Image) -> BaseCatalog:
     """
     Generate a photometric calibration catalog for SUMMER images

mirar/pipelines/winter/__init__.py

@@ -1 +1,4 @@
+"""
+Pipeline for WINTER data
+"""
 from mirar.pipelines.winter.winter_pipeline import WINTERPipeline