better astrometry (#454)

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

mirar/catalog/vizier/base_vizier_catalog.py

@@ -11,6 +11,14 @@
 from astroquery.vizier import Vizier
 
 from mirar.catalog.base_catalog import BaseCatalog
+from mirar.errors import ProcessorError
+
+
+class VizierError(ProcessorError):
+    """
+    Class for errors in Vizier catalog
+    """
+
 
 logger = logging.getLogger(__name__)
 
@@ -96,6 +104,14 @@ def get_catalog(self, ra_deg: float, dec_deg: float) -> astropy.table.Table:
 
         else:
             table = query[0]
+            logger.debug(f"Table columns are: {table.colnames}")
+            if self.get_mag_key() not in table.colnames:
+                err = (
+                    f"Magnitude column {self.get_mag_key()} not found in table."
+                    f"Available options are : {table.colnames}"
+                )
+                raise VizierError(err)
+
             table["ra"] = table[self.ra_key]
             table["dec"] = table[self.dec_key]
             table["magnitude"] = table[self.get_mag_key()]

mirar/pipelines/winter/blocks.py

@@ -2,7 +2,7 @@
 Module for WINTER data reduction
 """
 from mirar.io import open_fits
-from mirar.paths import BASE_NAME_KEY, FITS_MASK_KEY
+from mirar.paths import FITS_MASK_KEY
 from mirar.pipelines.winter.config import (
     sextractor_anet_config,
     sextractor_autoastrometry_config,
@@ -12,6 +12,7 @@
 from mirar.pipelines.winter.generator import (
     scamp_config_path,
     winter_astrometric_catalog_generator,
+    winter_astrostat_catalog_purifier,
     winter_photometric_catalog_generator,
     winter_reference_generator,
 )
@@ -52,7 +53,7 @@
     ImageSaver(output_dir_name="stacked_ref"),
 ]
 
-BOARD_ID = 4
+BOARD_ID = 2
 TARGET_NAME = "m39"
 split = [
     MultiExtParser(
@@ -229,6 +230,7 @@
         search_radius_deg=1.0,
         sextractor_config_path=sextractor_autoastrometry_config["config_path"],
         use_weight=True,
+        timeout=30,
     ),
     ImageSaver(output_dir_name=f"anet_{BOARD_ID}"),
     Sextractor(
@@ -261,7 +263,7 @@
     ImageSaver(output_dir_name="pre_anet"),
     AstrometryNet(
         output_sub_dir="anet",
-        scale_bounds=[10, 20],
+        scale_bounds=[15, 23],
         scale_units="amw",
         use_sextractor=True,
         parity="neg",
@@ -275,18 +277,31 @@
         write_regions_bool=True,
         output_sub_dir="scamp",
     ),
+    AstrometryStatsWriter(
+        ref_catalog_generator=winter_photometric_catalog_generator,
+        image_catalog_purifier=winter_astrostat_catalog_purifier,
+        write_regions=True,
+        cache=True,
+        crossmatch_radius_arcsec=5.0,
+    ),
+    ImageSaver(output_dir_name="anet"),
+    # Sextractor(
+    #     **sextractor_autoastrometry_config,
+    #     write_regions_bool=True,
+    #     output_sub_dir="scamp",
+    # ),
     ImageDebatcher(),
     ImageBatcher(["BOARD_ID", "FILTER", "TARGNAME", "SUBCOORD"]),
-    Scamp(
-        temp_output_sub_dir="scamp",
-        ref_catalog_generator=winter_astrometric_catalog_generator,
-        scamp_config_path=scamp_config_path,
-        cache=False,
-    ),
+    # Scamp(
+    #     temp_output_sub_dir="scamp",
+    #     ref_catalog_generator=winter_astrometric_catalog_generator,
+    #     scamp_config_path=scamp_config_path,
+    #     cache=False,
+    # ),
     Swarp(
         swarp_config_path=swarp_config_path,
         calculate_dims_in_swarp=True,
-        include_scamp=True,
+        include_scamp=False,
         subtract_bkg=False,
         cache=False,
         center_type="ALL",
@@ -478,8 +493,15 @@
 full_commissioning_all_boards = (
     load_unpacked + dark_cal_all_boards + flat_cal_all_boards + process_proc_all_boards
 )
-# commissioning_split = load_all_boards + split_images + process + \
-# process_proc_all_boards + photcal
+commissioning_split_single_board = (
+    log
+    + split_images
+    + dark_cal_all_boards
+    + flat_cal_all_boards
+    + process_proc_all_boards
+    + photcal
+)
+
 commissioning_split = (
     load_unpacked
     + split_images

mirar/pipelines/winter/config/files/astrom_anet.sex

@@ -5,18 +5,18 @@
 #-------------------------------- Catalog ------------------------------------
 
 CATALOG_NAME     test.cat       # name of the output catalog
-CATALOG_TYPE     ASCII_HEAD     # NONE,ASCII,ASCII_HEAD, ASCII_SKYCAT,
+CATALOG_TYPE     FITS_1.0       # NONE,ASCII,ASCII_HEAD, ASCII_SKYCAT,
                                 # ASCII_VOTABLE, FITS_1.0 or FITS_LDAC
-PARAMETERS_NAME  /Users/robertstein/Code/mirar/mirar/pipelines/winter/config/files/astrom.param  # name of the file containing catalog contents
+PARAMETERS_NAME  /Users/robertstein/Code/mirar/mirar/pipelines/winter/config/files/astrom_anet.param  # name of the file containing catalog contents
 
 #------------------------------- Extraction ----------------------------------
 
 DETECT_TYPE      CCD            # CCD (linear) or PHOTO (with gamma correction)
 DETECT_MINAREA   5              # minimum number of pixels above threshold
 THRESH_TYPE	 RELATIVE	# threshold type: RELATIVE (in sigmas) or ABSOLUTE (in ADUs)
 
-DETECT_THRESH    3            # <sigmas> or <threshold>,<ZP> in mag.arcsec-2
-ANALYSIS_THRESH  3            # <sigmas> or <threshold>,<ZP> in mag.arcsec-2
+DETECT_THRESH    1.5            # <sigmas> or <threshold>,<ZP> in mag.arcsec-2
+ANALYSIS_THRESH  1.5            # <sigmas> or <threshold>,<ZP> in mag.arcsec-2
 
 FILTER           Y              # apply filter for detection (Y or N)?
 FILTER_NAME      /Users/robertstein/Code/mirar/mirar/pipelines/winter/config/files/default.conv   # name of the file containing the filter
@@ -68,13 +68,13 @@ PIXEL_SCALE      0           # size of pixel in arcsec (0=use FITS WCS info)
 
 #------------------------- Star/Galaxy Separation ----------------------------
 
-SEEING_FWHM      2            # stellar FWHM in arcsec
+SEEING_FWHM      1.2            # stellar FWHM in arcsec
 STARNNW_NAME     /Users/robertstein/Code/mirar/mirar/pipelines/winter/config/files/default.nnw    # Neural-Network_Weight table filename
 
 #------------------------------ Background -----------------------------------
 
-BACK_SIZE        256             # Background mesh: <size> or <width>,<height>
-BACK_FILTERSIZE  6              # Background filter: <size> or <width>,<height>
+BACK_SIZE        64             # Background mesh: <size> or <width>,<height>
+BACK_FILTERSIZE  3              # Background filter: <size> or <width>,<height>
 
 BACK_TYPE        AUTO           # AUTO or MANUAL
 BACKPHOTO_TYPE   LOCAL        # can be GLOBAL or LOCAL

mirar/pipelines/winter/generator.py

@@ -16,6 +16,9 @@
 from mirar.pipelines.wirc.wirc_files import sextractor_astrometry_config
 from mirar.processors.astromatic.sextractor.sextractor import Sextractor
 from mirar.processors.astromatic.swarp.swarp import Swarp
+from mirar.processors.base_catalog_xmatch_processor import (
+    default_image_sextractor_catalog_purifier,
+)
 from mirar.processors.photcal import PhotCalibrator
 from mirar.processors.sqldatabase.base_model import BaseDB
 from mirar.references.local import RefFromPath
@@ -93,6 +96,16 @@ def winter_reference_image_resampler(**kwargs) -> Swarp:
     )
 
 
+def winter_astrostat_catalog_purifier(catalog: Table, image: Image) -> Table:
+    """
+    Default function to purify the photometric image catalog
+    """
+
+    return default_image_sextractor_catalog_purifier(
+        catalog, image, edge_width_pixels=0, fwhm_threshold_arcsec=20.0
+    )
+
+
 def winter_photometric_catalog_generator(image: Image) -> Gaia2Mass | PS1:
     """
     Function to crossmatch WIRC to GAIA/2mass for photometry
@@ -102,7 +115,9 @@ def winter_photometric_catalog_generator(image: Image) -> Gaia2Mass | PS1:
     """
     filter_name = image["FILTER"]
     search_radius_arcmin = (
-        np.max([image["NAXIS1"], image["NAXIS2"]]) * np.abs(image["CD1_1"]) * 60
+        np.max([image["NAXIS1"], image["NAXIS2"]])
+        * np.max([np.abs(image["CD1_1"]), np.abs(image["CD1_2"])])
+        * 60
     ) / 2.0
 
     if filter_name in ["J", "H"]:
@@ -119,31 +134,18 @@ def winter_photometric_catalog_generator(image: Image) -> Gaia2Mass | PS1:
             min_mag=10,
             max_mag=20,
             search_radius_arcmin=search_radius_arcmin,
-            filter_name=filter_name,
+            filter_name=filter_name.lower(),
         )
 
 
 def winter_ref_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]
+    return default_image_sextractor_catalog_purifier(
+        catalog, image, edge_width_pixels=100, fwhm_threshold_arcsec=4.0
+    )
 
 
 def winter_reference_phot_calibrator(image: Image, **kwargs) -> PhotCalibrator:

mirar/pipelines/winter/load_winter_image.py

@@ -49,7 +49,7 @@ def load_raw_winter_image(path: str | Path) -> tuple[np.array, astropy.io.fits.H
     logger.debug(f"Loading {path}")
     with fits.open(path) as img:
         # pylint: disable=E1101
-        data = img[0].data
+        data = img[0].data  # pylint: disable=no-member
         header = img[0].header
 
         header["UNIQTYPE"] = f"{header['OBSTYPE']}_{header['BOARD_ID']}"
@@ -136,7 +136,7 @@ def load_raw_winter_image(path: str | Path) -> tuple[np.array, astropy.io.fits.H
             pass
 
         if header["BOARD_ID"] == 2:
-            data[1085:, :] = np.nan
+            data[1060:, :] = np.nan
             data[:, 1970:] = np.nan
             data[:55, :] = np.nan
             data[:, :20] = np.nan
@@ -181,8 +181,8 @@ def load_proc_winter_image(path: str | Path) -> tuple[np.array, astropy.io.fits.
     """
     logger.debug(f"Loading {path}")
     with fits.open(path) as img:
-        data = img[0].data
-        header = img[0].header
+        data = img[0].data  # pylint: disable=no-member
+        header = img[0].header  # pylint: disable=no-member
         if "weight" in path:
             header["OBSTYPE"] = "weight"
 
@@ -199,8 +199,8 @@ def load_stacked_winter_image(
     """
     logger.info(f"Loading {path}")
     with fits.open(path) as img:
-        data = img[0].data
-        header = img[0].header
+        data = img[0].data  # pylint: disable=no-member
+        header = img[0].header  # pylint: disable=no-member
         if "weight" in path:
             header["OBSTYPE"] = "weight"
 

mirar/pipelines/winter/winter_pipeline.py

@@ -16,6 +16,7 @@
     commissioning_proc,
     commissioning_reduce,
     commissioning_split,
+    commissioning_split_single_board,
     commissioning_stack,
     full_commissioning,
     full_commissioning_all_boards,
@@ -54,6 +55,7 @@ class WINTERPipeline(Pipeline):
         "commissioning_split": commissioning_split,
         "unpack_subset": unpack_subset,
         "unpack_all": unpack_all,
+        "commissioning_split_single_board": commissioning_split_single_board,
     }
 
     gain = 1.0

mirar/pipelines/winter/write_obslog.py

@@ -0,0 +1,49 @@
+"""
+Module to write an observation log for a directory of images
+"""
+import argparse
+from glob import glob
+
+from astropy.io import fits
+from astropy.table import Table
+
+from mirar.paths import get_output_dir
+
+
+def write_observation_log(image_dir):
+    """
+    Write an observation log for a directory of images
+    """
+    filelist = glob(image_dir.as_posix() + "/*.fits")
+    keywords = [
+        "FILTER",
+        "EXPTIME",
+        "OBSTYPE",
+        "TARGNAME",
+        "RADEG",
+        "DECDEG",
+        "BOARD_ID",
+    ]
+
+    all_entries = []
+    for filename in filelist:
+        header = fits.getheader(filename)
+        entries = [header[key] for key in keywords if key in header]
+        found_keys = [key for key in keywords if key in header]
+        all_entries.append(entries)
+
+    tab = Table(rows=all_entries, names=found_keys)
+    tab["filename"] = [filename.split("/")[-1] for filename in filelist]
+    tab.write(f"{directory}/obslog.csv", format="ascii.csv", overwrite=True)
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description="Fix headers of images")
+    parser.add_argument("night", type=str, help="Night of observation")
+    parser.add_argument(
+        "--sub_dir", type=str, default="raw", help="Subdirectory of images"
+    )
+    args = parser.parse_args()
+
+    directory = get_output_dir(args.sub_dir, f"winter/{args.night}")
+    write_observation_log(directory)