From 61345f29453bf848e420741a57786b767d08c124 Mon Sep 17 00:00:00 2001 From: Saarah Hall <66533667+saarahha@users.noreply.github.com> Date: Mon, 26 Jun 2023 16:31:41 -0500 Subject: [PATCH] Better photometry for SEDMv2 (and beyond...) (#426) --- mirar/pipelines/sedmv2/blocks.py | 13 ++++++-- mirar/processors/photcal.py | 37 +++++++++++++++++++--- mirar/processors/utils/multi_ext_parser.py | 9 +++++- 3 files changed, 50 insertions(+), 9 deletions(-) diff --git a/mirar/pipelines/sedmv2/blocks.py b/mirar/pipelines/sedmv2/blocks.py index eb9287a6e..6cd106120 100644 --- a/mirar/pipelines/sedmv2/blocks.py +++ b/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), diff --git a/mirar/processors/photcal.py b/mirar/processors/photcal.py index f9642c624..ac40fde75 100644 --- a/mirar/processors/photcal.py +++ b/mirar/processors/photcal.py @@ -107,7 +107,9 @@ 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: @@ -115,14 +117,23 @@ def calculate_zeropoint(self, ref_cat_path: Path, img_cat_path: Path) -> list[di 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 diff --git a/mirar/processors/utils/multi_ext_parser.py b/mirar/processors/utils/multi_ext_parser.py index 643587d51..b3dc7dd51 100644 --- a/mirar/processors/utils/multi_ext_parser.py +++ b/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