midwayfalts

midwayfalts

Split up WINTER generator

Alphabetise error stack

More sorting

Select dome flats

Select dome flats

working

Nanfill

midwayfalts

midwayfalts

Split up WINTER generator

More sorting

Select dome flats

Select dome flats

working

Nanfill

mirar/pipelines/winter/blocks.py

@@ -39,6 +39,8 @@
 from mirar.pipelines.winter.constants import NXSPLIT, NYSPLIT
 from mirar.pipelines.winter.generator import (
     mask_stamps_around_bright_stars,
+    select_winter_dome_flats_images,
+    select_winter_flat_images,
     select_winter_sky_flat_images,
     winter_anet_sextractor_config_path_generator,
     winter_astrometric_ref_catalog_generator,
@@ -143,6 +145,7 @@
     ImageSaver,
     ImageSelector,
     MEFLoader,
+    NanFiller,
 )
 from mirar.processors.utils.cal_hunter import CalHunter
 from mirar.processors.xmatch import XMatch
@@ -304,6 +307,19 @@
 
 load_unpacked = [
     ImageLoader(input_sub_dir="raw_unpacked", input_img_dir=base_output_dir),
+    ImageSelector(
+        (
+            "TARGET",
+            [
+                "flat",
+                "dark",
+                "request_2024A002_2024_05_01_17_29_20.db_3",
+                "domeflat",
+                "FLAT",
+            ],
+        )
+    ),
+    # FIXME: remove the above line and uncomment the below line
     ImageRebatcher("UTCTIME"),
     CSVLog(
         export_keys=[
@@ -322,6 +338,7 @@
             "T_ROIC",
             "FIELDID",
             "MEDCOUNT",
+            "DITHGRP",
         ]
     ),
 ]
@@ -350,15 +367,34 @@
 ]
 
 flat_calibrate = [
-    ImageRebatcher(["BOARD_ID", "FILTER"]),
-    ImageSelector((OBSCLASS_KEY, ["science"])),
+    ImageRebatcher(["SUBCOORD", "FILTER"]),
     FlatCalibrator(
         cache_sub_dir="calibration_flats",
-        select_flat_images=select_winter_sky_flat_images,
+        select_flat_images=select_winter_dome_flats_images,
         cache_image_name_header_keys=["FILTER", "BOARD_ID"],
+        flat_mode="pixel",
     ),
-    ImageSaver(output_dir_name="skyflatcal"),
-    ImageRebatcher(["BOARD_ID", "UTCTIME", "SUBCOORD"]),
+    ImageSelector((OBSCLASS_KEY, ["science"])),
+    ImageSaver(output_dir_name="domeflatcal"),
+    # FIXME Add in the old sky flat thing here post dome flat
+    ImageRebatcher(
+        [
+            "BOARD_ID",
+            "FILTER",
+            "SUBCOORD",
+            "GAINCOLT",
+            "GAINCOLB",
+            "GAINROW",
+            TARGET_KEY,
+        ]
+    ),
+    FlatCalibrator(
+        cache_sub_dir="sky_dither_flats",
+        select_flat_images=select_winter_sky_flat_images,
+        flat_mode="structure",
+    ),
+    ImageSaver(output_dir_name="allskyflatcal"),
+    ImageRebatcher([BASE_NAME_KEY]),
     Sextractor(
         **sextractor_astrometry_config,
         write_regions_bool=True,
@@ -379,7 +415,7 @@
 ]
 
 astrometry = [
-    ImageRebatcher(["UTCTIME", "BOARD_ID", "SUBCOORD"]),
+    ImageRebatcher([BASE_NAME_KEY]),
     AstrometryNet(
         output_sub_dir="anet",
         scale_bounds=[1.0, 1.3],
@@ -432,6 +468,9 @@
     ImageDebatcher(),
     CustomImageBatchModifier(winter_boardid_6_demasker),
     ImageBatcher(["BOARD_ID", "FILTER", "EXPTIME", TARGET_KEY, "SUBCOORD"]),
+    ImageSaver(output_dir_name="prestack"),
+    NanFiller(),
+    ImageSaver(output_dir_name="postnanfill"),
     Swarp(
         swarp_config_path=swarp_config_path,
         calculate_dims_in_swarp=True,

mirar/pipelines/winter/generator/__init__.py

@@ -25,6 +25,7 @@
 from mirar.pipelines.winter.generator.reduce import (
     mask_stamps_around_bright_stars,
     select_winter_dome_flats_images,
+    select_winter_flat_images,
     select_winter_sky_flat_images,
     winter_anet_sextractor_config_path_generator,
     winter_boardid_6_demasker,

mirar/pipelines/winter/generator/reduce.py

@@ -20,7 +20,7 @@
     OBSCLASS_KEY,
     REF_CAT_PATH_KEY,
     SATURATE_KEY,
-    get_output_dir,
+    get_output_dir
 )
 from mirar.pipelines.winter.config import sextractor_anet_config
 from mirar.pipelines.winter.fourier_bkg_model import subtract_fourier_background_model

mirar/pipelines/winter/load_winter_image.py

@@ -543,7 +543,7 @@ def get_raw_winter_mask(image: Image) -> np.ndarray:
         mask[1060:, :] = 1.0
 
         # Mask a vertical strip
-        mask[:, 998:1002] = 1.0
+        mask[:, 992:1002] = 1.0
 
         # Mask another vertical strip
         mask[:, 1266:1273] = 1.0

mirar/processors/flat.py

@@ -9,6 +9,15 @@
 from copy import copy
 
 import numpy as np
+<<<<<<< HEAD
+<<<<<<< HEAD
+from astropy.convolution import Tophat2DKernel, convolve_fft
+=======
+from astropy.convolution import Tophat2DKernel, convolve, convolve_fft
+>>>>>>> 60685c68 (working)
+=======
+from astropy.convolution import Tophat2DKernel, convolve_fft
+>>>>>>> 822f21d8 (Nanfill)
 
 from mirar.data import Image, ImageBatch
 from mirar.errors import ImageNotFoundError
@@ -26,6 +35,28 @@
 logger = logging.getLogger(__name__)
 
 
+def get_convolution(data: np.ndarray, kernel_width: int) -> np.ndarray:
+    """
+    Convolve data with a tophat kernel
+
+    :param data: Image data
+    :param kernel_width: Width of the kernel (pixels)
+    :return: Smoothed image
+    """
+    pad_top = np.array([data[0] for _ in range(kernel_width)])
+    pad_bottom = np.array([data[-1] for _ in range(kernel_width)])
+    extended = np.vstack([pad_top, data, pad_bottom])
+    pad_left = np.array([extended.T[0] for _ in range(kernel_width)])
+    pad_right = np.array([extended.T[-1] for _ in range(kernel_width)])
+    extended = np.hstack([pad_left.T, extended, pad_right.T])
+
+    tophat_kernel = Tophat2DKernel(kernel_width)
+    smooth_illumination = convolve_fft(
+        extended, tophat_kernel, nan_treatment="interpolate"
+    )[kernel_width:-kernel_width, kernel_width:-kernel_width]
+    return smooth_illumination
+
+
 class MissingFlatError(ImageNotFoundError):
     """
     Error for when a dark image is missing
@@ -61,6 +92,7 @@ def __init__(
         flat_nan_threshold: float = 0.0,
         select_flat_images: Callable[[ImageBatch], ImageBatch] = default_select_flat,
         flat_mask_key: str = None,
+        flat_mode: str = "median",
         **kwargs,
     ):
         super().__init__(*args, **kwargs)
@@ -71,6 +103,9 @@ def __init__(
         self.flat_nan_threshold = flat_nan_threshold
         self.select_cache_images = select_flat_images
         self.flat_mask_key = flat_mask_key
+        self.flat_mode = flat_mode
+        if not self.flat_mode in ["median", "pixel", "structure"]:
+            raise ValueError(f"Flat mode {self.flat_mode} not supported")
 
     def __str__(self) -> str:
         return "Creates a flat image, divides other images by this image."
@@ -152,6 +187,53 @@ def make_image(
 
         master_flat = np.nanmedian(flats, axis=2)
 
+        if self.flat_mode != "median":
+
+            if self.flat_mode == "pixel":
+                flatdata_norm_smooth = get_convolution(master_flat, 40)
+
+                pixel_variation = master_flat / flatdata_norm_smooth
+
+                # Clip outliers (they'll get worked out in stacking)
+                std = np.nanstd(pixel_variation)
+                sig = abs(pixel_variation - np.nanmedian(pixel_variation)) / std
+
+                mask = sig > 2.0
+
+                logger.debug(
+                    f"Masking {np.sum(mask)} pixels "
+                    f"out of {len(mask.flatten()) }in flat"
+                )
+                pixel_variation[mask] = np.nan
+                master_flat = pixel_variation / np.nanmedian(pixel_variation)
+
+            elif self.flat_mode == "structure":
+                flatdata_norm_smooth = get_convolution(master_flat, 100)
+                flatdata_norm_smooth[np.isnan(master_flat)] = np.nan
+
+                pixel_variation = master_flat / flatdata_norm_smooth
+
+                # Clip outliers (they'll get worked out in stacking)
+                std = np.nanstd(pixel_variation)
+                sig = abs(pixel_variation - np.nanmedian(pixel_variation)) / std
+
+                mask = sig > 3.0
+
+                logger.debug(
+                    f"Masking {np.sum(mask)} pixels "
+                    f"out of {len(mask.flatten()) }in flat"
+                )
+
+                pixel_variation[sig > 3.0] = np.nan
+                master_flat = (
+                    pixel_variation
+                    * flatdata_norm_smooth
+                    / np.nanmedian(pixel_variation)
+                )
+
+            else:
+                raise ValueError(f"Flat mode {self.flat_mode} not supported")
+
         master_flat_image = Image(master_flat, header=copy(images[0].get_header()))
         master_flat_image[COADD_KEY] = n_frames
 

mirar/processors/utils/__init__.py

@@ -17,3 +17,4 @@
     select_from_images,
 )
 from mirar.processors.utils.multi_ext_parser import MultiExtParser
+from mirar.processors.utils.nan_filler import NanFiller

mirar/processors/utils/nan_filler.py

@@ -0,0 +1,35 @@
+"""
+Module for filling NaN values in an image
+"""
+
+import logging
+
+import numpy as np
+
+from mirar.data import ImageBatch
+from mirar.processors.base_processor import BaseImageProcessor
+
+logger = logging.getLogger(__name__)
+
+
+class NanFiller(BaseImageProcessor):
+    """
+    Processor to fill image nans
+    """
+
+    base_key = "nanfill"
+
+    def __str__(self):
+        return "Processor to pad masked image regions with the median."
+
+    def _apply_to_images(
+        self,
+        batch: ImageBatch,
+    ) -> ImageBatch:
+
+        for image in batch:
+            data = image.get_data()
+            data[np.isnan(data)] = np.nanmedian(data)
+            image.set_data(data)
+
+        return batch