CR not removed by median stacking

By default, MEGARA DRP recipes use median combination of three or more images to produce an average exposure and eliminate cosmic rays. However, when exposure times are long, the same pixel can be hit by cosmic rays in multiple consecutive exposures, causing the contamination to persist in the median-combined image.

The most straightforward solution is to visually inspect the resulting image (e.g., the final_rss.fits file) and manually interpolate the affected pixels using neighboring data. For MEGARA RSS images, this approach is only reliable when performing 1D interpolation along the wavelength direction, since the spectra of neighboring fibers in this file do not correspond to neighboring spaxels in the field of view. Interactive programs such as tea-cleanest can facilitate this interpolation task.

However, removing residual cosmic ray pixels in the RSS image produce worse results than correcting the contamination before combining the signal from each fiber at a given wavelength. A more efficient approach is to use a procedure that detects and corrects cosmic ray-affected pixels across multiple exposures prior to extracting the resulting spectrum from each fiber.

Warning

The method described here is available in the numina crmasks. subpackage. We strongly recommend that users review that documentation thoroughly before applying it to MEGARA data. This method can be easily activated by modifying the default image combination settings in the MEGARA DRP.

Requirements and methodology

This method requires three or more equivalent exposures. It identifies pixels where cosmic rays have affected the signal in more than half of the available exposures—cases where median combination alone cannot remove the contamination.

Note

This method is likely to be useful primarily when reducing science images with long exposure times. In general, it is not recommended for the reduction of calibration images.

Overall procedure

The main workflow for applying the method is as follows:

  1. Prepare the observation result file: Create the YAML file that lists the individual exposures and specifies additional requirements for the reduction recipe associated with the reduction recipe MegaraCrDetection.

  2. Run the MEGARA DRP to generate the CR masks: The execution of the MEGARA DRP with this YAML file is interactive (see below), and the result is a FITS file named crmasks.fits, which contains several masks that can later be used to remove the cosmic rays using different strategies. This file must be copied to the data/ subdirectory.

  3. Run the MEGARA DRP to generate the reduced scientific result: Once the file crmasks.fits is placed in the data/ subdirectory, the MEGARA DRP can be run again to generate the reduced scientific result, using for that purpose the corresponding YAML file. During this step, the user must chose one of the available strategies for combining exposures and removing cosmic rays. Instead of using the default median combination, one must specify one of the following methods in the requirements section: mediancr, meancrt or meancr.

Note

The available combination methods, mediancr, meancrt, meancr and meancr2, produce different results:

  • mediancr: performs a median combination, replacing pixels suspected of being affected by cosmic rays in more than one exposure by the “minimum value” at each pixel across the available exposures.

  • meancrt: generates the mean combination using of a single mask that stores all the cosmic rays detected in all the individual exposures, replacing each masked pixels by the value obtained when using the mediancr method.

  • meancr: also performs a mean combination, but uses the individual cosmic ray masks specifically computed for each available exposure.

  • meancr2: makes use of meancr and interpolates residual cosmic ray pixels detected using the selected method.

Since the mean has a lower standard deviation than the median, the meancrt and meancr methods are generally preferable. Among these two, tests suggest that meancr tends to yield better results. However, it is recommended to try all three methods and compare the outputs to determine which works best for your specific dataset.

Preparing the observation result file

Let’s assume that we have three initially equivalent exposures. The process begins with the creation of a YAML file, e.g. 8_generate_crmasks.yaml, which lists the individual exposures to be used in the cosmic ray detection step.

 1id: 8_LR-R_crmasks_mm_pycosmic2
 2mode: MegaraCrDetection
 3instrument: MEGARA
 4frames:
 5 - 0003978527-20231115-MEGARA-MegaraLcbImage.fits
 6 - 0003978528-20231115-MEGARA-MegaraLcbImage.fits
 7 - 0003978529-20231115-MEGARA-MegaraLcbImage.fits
 8requirements:
 9  # Detection methods:
10  # mm_lacosmic | mm_pycosmic | mm_deepcr | mm_conn |
11  # lacosmic | pycosmic | deepcr | conn
12  crmethod: mm_pycosmic             # any of the above methods
13  use_auxmedian: False              # use aux. Median insteand of minimum value
14  interactive: True                 # display and pause between plots
15  flux_factor: none                 # none | auto | [list of values]
16  flux_factor_regions:              # list of regions, FITS criterium
17  # - [xmin, xmax, ymin, ymax]      # <-- format (no need to uncomment this)
18  apply_flux_factor_to: simulated   # original | simulated (data)
19  dilation: 0                       # global dilation (integer)
20  regions_to_be_skipped:            # list of regions, FITS criterium
21  # - [xmin, xmax, ymin, ymax]      # <-- format (no need to uncomment this)
22  pixels_to_be_flagged_as_cr:       # FITS criterium, first pixel is [1, 1]
23  # - [X, Y]                        # <-- format (no need to uncomment this)
24  pixels_to_be_ignored_as_cr:       # FITS criterium, first pixel is [1, 1]
25  # - [X, Y]                        # <-- format (no need to uncomment this)
26  pixels_to_be_replaced_by_local_median:  # FITS criterium
27  # - [X, Y, X_width, Y_width]      # <-- format (no need to uncomment this)
28  verify_cr: False                  # ask for confirmation
29  semiwindow: 15                    # to display residual CR hits
30  color_scale: zscale               # minmax | zscale
31  maxplots: -1                      # -1: all CR detected are plotted
32  #-----------------------------------------------------------------------------
33  # Parameters for L.A. Cosmic (van Dokkum 2001), as implemented in ccdproc
34  # (based in https://github.com/astropy/astroscrappy).
35  # https://ccdproc.readthedocs.io/en/latest/api/ccdproc.cosmicray_lacosmic.html
36  la_gain_apply: True    # Return gain-corrected data
37  la_sigclip: [5.0,3.0]  # Laplacian-to-noise limit for cosmic ray detection
38  la_sigfrac: [0.3,0.3]  # Fractional detection limit for neighboring pixels
39  la_objlim: [5.0,5.0]   # Contrast between Laplacian and fine struct. images
40  la_satlevel: 65535.0   # Saturation level of the image (ADU)
41  la_niter: 4            # Number of iterations to perform
42  la_sepmed: True        # Use the separable median filter (not the full filter)
43  la_cleantype: meanmask # median | medmask | meanmask | idw
44  la_fsmode: convolve    # median | convolve
45  la_psfmodel: gaussxy   # gauss | moffat | gaussx | gaussy | gaussxy
46  la_psffwhm_x: 2.5      # FWHM of the PSF (X direction) to generate the kernel
47  la_psffwhm_y: 2.5      # FWHM of the PSF (Y direction) to generate the kernel
48  la_psfsize: 11         # size of the kernel
49  la_psfbeta: 4.765      # Moffat beta parameter
50  la_verbose: True       # Print to the screen or not
51  la_padwidth: 10        # Padding to mitigate edge effects [not in original]
52  #-----------------------------------------------------------------------------
53  # Parameters for PyCosmic (Husemann et al. 2012)
54  # Note: gain, rdnoise and bias will be set to the general values
55  pc_sigma_det: [5.0,3.0] # Detection limit above the noise
56  pc_rlim: [1.2,1.2]      # Detection threshold
57  pc_iterations: 5        # Number of iterations
58  pc_fwhm_gauss_x: 2.5    # FWHM of the Gaussian smoothing kernel (X direction)
59  pc_fwhm_gauss_y: 2.5    # FWHM of the Gaussian smoothing kernel (Y direction)
60  pc_replace_box_x: 5     # Median box size (X axis) to estimate replacement
61  pc_replace_box_y: 5     # Median box size (Y axis) to estimate replacement
62  pc_replace_error: 1.0E6 # Error value for bad pixels
63  pc_increase_radius: 0   # Dilation (pixels)
64  pc_verbose: True        # Provide information during the processing
65  #-----------------------------------------------------------------------------
66  # Parameters for deepCR (Zhang & Bloom 2020)
67  dc_mask: ACS-WFC        # ACS-WFC | WFC3-UVIS
68  dc_threshold: 0.5       # Probability threshold to be applied
69  dc_verbose: True        # Provide information during the processing
70  #-----------------------------------------------------------------------------
71  # Parameters for Cosmic-CoNN (Xu et al. 2023)
72  nn_model: ground_imaging # ground_imaging | NRES | HST_ACS_WFC
73  nn_threshold: 0.5        # Probability threshold to be applied
74  nn_verbose: True         # Provide information during the processing
75  #-----------------------------------------------------------------------------
76  # Parameters for the Median-Minimum method (Cardiel et al. 2026)
77  mm_cr_coincidences: 2      # number of CR coincidences to be sought
78  mm_photon_distribution: poisson  # poisson | nbinom
79  mm_nbinom_shape: 1000      # shape parameter for negative binomial (nbinom)
80  mm_synthetic: single       # median | single (to simulate M.M. diagram)
81  mm_hist2d_min_neighbors: 8 # minimum number of neighbors in hist2d bins
82  mm_ydiag_max: 0            # maximum Y value in hist2d (0=auto)
83  mm_dilation: 1             # dilation around suspected pixels (integer)
84  # rectangular region to determine offsets between individual images:
85  mm_xy_offsets:             # XYoffsets between individual images (pixels)
86  # - [xoffset, yoffset]     # pair of values for each individual exposure
87  mm_crosscorr_region: null  # [xmin, xmax, ymin, ymax] FITS criterium | null
88  mm_boundary_fit: spline    # piecewise | spline
89  mm_knots_splfit: 5         # used only for spline fit
90  mm_fixed_points_in_boundary:
91  # - [x, y, weight]  # weights are optional(ignored) for spline(piecewise) fit
92  mm_nsimulations: 10                 # simulations for each set of images
93  mm_niter_boundary_extension: 10     # iterations to extend the spline fit
94  mm_weight_boundary_extension: 2.0   # multiplicative increment in weigth
95  mm_threshold_rnoise: 3.0            # threshold value (in readout units)
96  mm_minimum_max2d_rnoise: 5.0        # minimum max2d (in readout units)
97  mm_seed: 123456                     # random seed for reproducibility
98  #----------------------------------------------------------------------------

Note that in this file is important to specify:

  • mode: MegaraCrDetection, so the MEGARA DRP will know that the 8_generate_crmasks.yaml file is intended for generating the CR masks.

  • frames: one should provide here the list of the individual equivalent exposures.

  • requirements: this section contains a large set of parameters whose values define the strategy used for the cosmic-ray detection. This section is identical to the requirements section of the configuration YAML file employed by numina-crmasks, which is described in detail in this link.

Note

In addition to the extensive list of requirements displayed above, inherited from the requirements section of the input YAML file used by numina-crmasks, this recipe also accepts two additional requirements useful for debugging purposes:

  • save_preprocessed: True (False by default): if True, the following images are saved in the corresponding work subdirectory:

    • preprocessed_1.fits, preprocessed_2.fits and preprocessed_3.fits: the preprocessed individual exposures after overscan correction, trimming, bias subtraction, gain correction and flatfielding (the images data values are in electrons).

    • preprocessed_median.fits: the median combination of the thre previous preprocessed individual images.

  • save_postprocessed: True (False by default): if True, the following images are also saved in the corresponding work subdirectory: postprocessed_mean.fits, postprocessed_min.fits, postprocessed_median.fits, postprocessed_mediancr.fits, postprocessed_meancrt.fits, postprocessed_meancr.fits, and postprocessed_meancr2.fits.

    These files correspond to the different combinations allowed by numina-crmasks, as explained in this link. Note that the contents of postprocessed_median.fits should be identical to preprocessed_median.fits.

Computing the CR masks

Coincident cosmic rays in the median combination

The recipe is run by doing:

(megara) $ numina run 8_generate_crmasks.yaml -r control.yaml

From this point, follow the description given in these examples, which helps understanding the impact of using different parameters in the cosmic-ray pixel detection process. Important: in those examples the generation of the CR masks is carried out by executing the command-line script numina-crmasks, which makes use of a slightly different input YAML file. However, the code executed is the same as in the MegaraCrDetection recipe, so the resulting output is identical.

Output file crmasks.fits

The output file crmasks.fits is generated in the results/ subdirectory, and contains several extensions:

(megara) $ cd obsid8_LcbImage_LR-R_crmasks_results

(megara) $ fitsinfo crmasks.fits

Filename: crmasks.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU     121   ()
  1  MEDIANCR      1 ImageHDU         8   (4096, 4112)   uint8
  2  MEANCRT       1 ImageHDU         8   (4096, 4112)   uint8
  3  CRMASK1       1 ImageHDU         8   (4096, 4112)   uint8
  4  CRMASK2       1 ImageHDU         8   (4096, 4112)   uint8
  5  CRMASK3       1 ImageHDU         8   (4096, 4112)   uint8
  6  MEANCR        1 ImageHDU         8   (4096, 4112)   uint8
  7  AUXCLEAN      1 ImageHDU         8   (4096, 4112)   float32

The primary HDU does not contain image data but includes keywords that document the parameters used to generate the masks.

The next extensions store the actual masks computed by the program:

  • MEDIANCR: Mask for pixels suspected of being affected by two cosmic rays in the median combination.

  • MEANCRT: Mask for pixels suspected of being affected by a cosmic ray in the mean combination.

  • CRMASK1, CRMASK2, CRMASK3: Masks for pixels suspected of being affected by a cosmic ray in each of the individual exposures.

  • MEANCR: Mask computed using the chosen CR detection method on the masked mean combination.

In all these masks, pixels suspected of being affected by a cosmic ray are set to 1, while the rest of the pixels are set to 0.

Finally, the extension AUXCLEAN contains the corrected image generated by the chosen CR detection algorithm.

Do not forget to copy the file crmasks.fits to the data/ subdirectory before running the MEGARA DRP to generate the reduced scientific result.

(megara) $ cp crmasks.fits ../data/

Applying the masks

Once the crmasks.fits file has been generated, it can be used to remove the suspected cosmic rays from the combined image. Several options are available for this step, depending on the value of method specified in the requirements section of the YAML file used to produce the reduced scientific result. Rather than assuming the combination method is median (the default in MEGARA DRP), the user should explicitly specify one of the following methods: mediancr, meancrt, meancr or meancr2.

Method mediancr

For example, when using the MegaraLcbImage recipe, the corresponding YAML file, 8_LcbImage.yaml should define the method and crmasks parameters in the requirements section, as follows:

 1id: 8_LcbImage_LR-R
 2mode: MegaraLcbImage
 3instrument: MEGARA
 4frames:
 5 - 0003978527-20231115-MEGARA-MegaraLcbImage.fits
 6 - 0003978528-20231115-MEGARA-MegaraLcbImage.fits
 7 - 0003978529-20231115-MEGARA-MegaraLcbImage.fits
 8requirements:
 9  method: mediancr
10  method_kwargs:
11    apply_flux_factor: False
12    use_auxmedian: False
13  crmasks: ../data/crmasks.fits
14  extraction_offset: [+3]
15  ignored_sky_bundles: [96] #[93, 94, 96, 97, 99, 100]

Note that the crmasks parameter must be explicitly specified, and it should point to the crmasks.fits file generated in the previous step. This requirement also gives users the flexibility to rename the file and test different versions of it (each potentially created with varying detection thresholds) to evaluate the impact of using different masks.

(megara) $ numina run 8_LcbImage.yaml -r control.yaml

When using method: mediancr, the MEGARA DRP reads the mask stored in the MEDIANCR extension of the crmasks.fits file. The masked pixels (those suspected of being affected by a cosmic ray in more than one exposure) are replaced with the corresponding min2d value at each pixel. As a result, the final image is identical to median2d, except for the corrected pixels.

Note that under the label method_kwargs: there are several parameters that control how the information from the crmasks.fits file is used:

  • apply_flux_factor: When set to True, the factors listed in the FLUXF1, FLUXF2, etc. keywords of the primary HDU in the crmasks.fits file are used to modify the signal of each individual exposure before combining them and removing cosmic rays. If this parameter is set to False, an effective factor of 1.0 is applied (i.e., no change is introduced), regardless of the values in the FITS header of the crmasks.fits file.

  • use_auxmedian:: When set to True, the masked pixels are replaced with the corrected value obtained using chosen CR detection algorithm, rather than using the “minimum value” in each pixel.

Method meancrt

 1id: 8_LcbImage_LR-R
 2mode: MegaraLcbImage
 3instrument: MEGARA
 4frames:
 5 - 0003978527-20231115-MEGARA-MegaraLcbImage.fits
 6 - 0003978528-20231115-MEGARA-MegaraLcbImage.fits
 7 - 0003978529-20231115-MEGARA-MegaraLcbImage.fits
 8requirements:
 9  method: meancrt
10  method_kwargs:
11    apply_flux_factor: False
12    use_auxmedian: False
13  crmasks: ../data/crmasks.fits
14  extraction_offset: [+3]
15  ignored_sky_bundles: [96] #[93, 94, 96, 97, 99, 100]

When using method: meancrt, the MEGARA DRP reads the mask stored in the MEANCRT extension of the file crmasks.fits. The masked pixels, which are suspected of being affected by a cosmic ray, are replaced by the corresponding values when using method: mediancr.

As a result, the final image is identical to mean2d, except for the corrected pixels, where the data is sourced from the image produced by the mediancr method.

The parameters under the method_kwargs: label play the same role as described above for the mediancr combination method.

Method meancr

 1id: 8_LcbImage_LR-R
 2mode: MegaraLcbImage
 3instrument: MEGARA
 4frames:
 5 - 0003978527-20231115-MEGARA-MegaraLcbImage.fits
 6 - 0003978528-20231115-MEGARA-MegaraLcbImage.fits
 7 - 0003978529-20231115-MEGARA-MegaraLcbImage.fits
 8requirements:
 9  method: meancr
10  method_kwargs:
11    apply_flux_factor: False
12    use_auxmedian: False
13  crmasks: ../data/crmasks.fits
14  extraction_offset: [+3]
15  ignored_sky_bundles: [96] #[93, 94, 96, 97, 99, 100]

When using method: meancr, the MEGARA DRP reads the masks stored in the CRMASK1, CRMASK2 and CRMASK3 extensions of the file crmasks.fits. The program uses NumPy masked arrays to construct a 3D stack of the individual exposures, where each image is represented as a masked array based on its corresponding mask. This allows the program to compute the mean along the axis corresponding to the image number, excluding the masked pixels.

In cases where a pixel is masked in all individual exposures, the corresponding pixel in the combined image is set to either the corresponding min2d value (if use_auxmedian: False) or to the value from the median combination computed using the mediancr method (if use_auxmedian: True).

The parameters under the method_kwargs: label play the same role as described above for the mediancr combination method.

Method meancr2

 1id: 8_LcbImage_LR-R
 2mode: MegaraLcbImage
 3instrument: MEGARA
 4frames:
 5 - 0003978527-20231115-MEGARA-MegaraLcbImage.fits
 6 - 0003978528-20231115-MEGARA-MegaraLcbImage.fits
 7 - 0003978529-20231115-MEGARA-MegaraLcbImage.fits
 8requirements:
 9  method: meancr
10  method_kwargs:
11    apply_flux_factor: False
12    use_auxmedian: False
13  crmasks: ../data/crmasks.fits
14  extraction_offset: [+3]
15  ignored_sky_bundles: [96] #[93, 94, 96, 97, 99, 100]

This is a refined version of the median combination. The method begins with the result from the previous case (method meancr), which is then cleaned using a mask generated by the same method chosen to detect cosmic rays in single exposures. Detected cosmic-ray pixels are replaced by either the corresponding min2d values (if use_auxmedian: False) or by values from the median combination computed using the mediancr method (if use_auxmedian: True).

The parameters under the method_kwargs: label play the same role as described above for the mediancr combination method.