Better photometry for SEDMv2 (and beyond...) (#426)

mirar/pipelines/sedmv2/blocks.py

@@ -39,6 +39,7 @@
 from mirar.processors.reference import ProcessReference
 from mirar.processors.utils import (
     ImageBatcher,
+    ImageDebatcher,
     ImageLoader,
     ImageSaver,
     ImageSelector,
@@ -53,7 +54,7 @@
 )
 
 load_raw = [
-    MultiExtParser(input_sub_dir="raw/mef/"),
+    MultiExtParser(input_sub_dir="raw/mef/", skip_first=True),
     ImageLoader(load_image=load_raw_sedmv2_image),
 ]
 
@@ -78,6 +79,7 @@
 ]  # pylint: disable=duplicate-code
 
 reduce = [
+    MaskPixels(mask_path=sedmv2_mask_path),
     BiasCalibrator(),
     ImageSelector(("OBSTYPE", ["FLAT", "SCIENCE"])),
     ImageBatcher(split_key="filter"),
@@ -92,7 +94,7 @@
         downsample=2,
         timeout=900,
     ),
-    MaskPixels(mask_path=sedmv2_mask_path),
+    ImageSaver(output_dir_name="masked", write_mask=True),
     Sextractor(
         output_sub_dir="sextractor",
         checkimage_name=None,
@@ -176,17 +178,22 @@
 parse_transient = [ImageSelector(("SOURCE", ["transient", "None"]))]
 
 resample_transient = [
+    ImageDebatcher(),
+    ImageBatcher(split_key="origname"),  # reaches for files coming from the same MEF
     Swarp(
+        cache=True,
         swarp_config_path=swarp_config_path,
         include_scamp=False,
         combine=True,
+        calculate_dims_in_swarp=True,
     ),
     ImageSaver(
         output_dir_name="resampled", write_mask=True
     ),  # pylint: disable=duplicate-code
 ]
 
-process_transient = parse_transient + reduce + resample_transient + calibrate
+# process_transient = parse_transient + reduce + resample_transient + calibrate
+process_transient = reduce + resample_transient + calibrate
 
 subtract = [
     ImageBatcher(split_key=BASE_NAME_KEY),

mirar/processors/photcal.py

@@ -107,22 +107,33 @@ def get_phot_output_dir(self):
         """
         return get_output_dir(self.temp_output_sub_dir, self.night_sub_dir)
 
-    def calculate_zeropoint(self, ref_cat_path: Path, img_cat_path: Path) -> list[dict]:
+    def calculate_zeropoint(
+        self, ref_cat_path: Path, img_cat_path: Path, img_filt
+    ) -> list[dict]:
         """
         Function to calculate zero point from two catalogs
         Args:
             ref_cat_path: Path to reference ldac catalog
             img_cat_path: Path to image ldac catalog
         Returns:
         """
-        ref_cat = get_table_from_ldac(ref_cat_path)
+        ref_cat_with_flagged = get_table_from_ldac(ref_cat_path)
         img_cat = get_table_from_ldac(img_cat_path)
 
-        if len(ref_cat) == 0:
+        if len(ref_cat_with_flagged) == 0:
             err = "No sources found in reference catalog"
             logger.error(err)
             raise PhotometryReferenceError(err)
 
+        if str(ref_cat_path).split(".")[-2] == "ps1":
+            # this reference catalog is from ps1
+            # remove sources with SATURATED flag
+            ref_cat = self.remove_sat_ps1(ref_cat_with_flagged, img_filt)
+
+        else:
+            # reference not ps1, no flags to check
+            ref_cat = ref_cat_with_flagged
+
         ref_coords = SkyCoord(ra=ref_cat["ra"], dec=ref_cat["dec"], unit=(u.deg, u.deg))
         clean_mask = (
             (img_cat["FLAGS"] == 0)
@@ -285,8 +296,9 @@ def _apply_to_images(
             image["FWHM_MED"] = fwhm_med
             image["FWHM_STD"] = fwhm_std
 
-            zp_dicts = self.calculate_zeropoint(ref_cat_path, temp_cat_path)
-
+            zp_dicts = self.calculate_zeropoint(
+                ref_cat_path, temp_cat_path, image.header["FILTER"]
+            )
             aperture_diameters = []
             zp_values = []
             for zpvals in zp_dicts:
@@ -440,3 +452,18 @@ def get_sextractor_apertures(self) -> list[float]:
         line = aperture_lines[0].replace("PHOT_APERTURES", " ").split("#")[0]
 
         return [float(x) for x in line.split(",") if x not in [""]]
+
+    def remove_sat_ps1(self, catalog, filt: str):
+        """
+        remove ps1 sources flagged as "SATURATED"
+        """
+        logger.info(f"original ps1 table length: {len(catalog)}")
+        logger.info("removing ps1 sources with SATURATED flag...")
+        sat_flag = 4096  # SATURATED value
+        column = catalog[str(filt) + "Flags"]
+        check = (column & sat_flag) / sat_flag
+        # check != 0 means this flag is there
+        # check == 0 means this flag is not there
+        clean_cat = catalog[np.where(check == 0)[0]]
+        logger.info(f"found {len(clean_cat)} columns without this flag \n")
+        return clean_cat

mirar/processors/utils/multi_ext_parser.py

@@ -32,11 +32,13 @@ def __init__(
         input_sub_dir: str = RAW_IMG_SUB_DIR,
         input_img_dir: str = base_raw_dir,
         load_image: Callable[[str], [np.ndarray, astropy.io.fits.Header]] = open_fits,
+        skip_first: bool = False,
     ):
         super().__init__()
         self.input_sub_dir = input_sub_dir
         self.input_img_dir = input_img_dir
         self.load_image = load_image
+        self.skip_first = skip_first
 
     def __str__(self):
         return (
@@ -63,7 +65,12 @@ def parse(self, path: str) -> list:
             # zip hdr0's values and comments
             zipped = list(zip(hdr0.values(), hdr0.comments))
             # combining main header (hdr0) with extension header
-            for ext in range(1, num_ext):
+            if self.skip_first:
+                start = 2
+                logger.info("Ignoring first extension frame")
+            else:
+                start = 1
+            for ext in range(start, num_ext):
                 data = hdu[ext].data
                 hdrext = hdu[ext].header