One way of producing flat field images is to use the science images, combined in a way that eliminates astronomical sources. This provides an exact match to the spectrum of the night sky, since the night sky is the source of light. However, the night sky is dark, so the counts in individual images is low. Many images must be combined to generate a flat with low noise.
from pathlib import Path
import numpy as np
from astropy.nddata import CCDData
from astropy import units as u
from astropy.stats import mad_std
from photutils import detect_threshold, detect_sources, source_properties
import ccdproc as ccdp
There are a few circumstances in which producing a sky flat is difficult or impossible:
roughly the same pixels in all of the images. In this case, there is no way to produce a good flat. If several fields of view are observed this should not be an issue.
of view. In this case there is likely to be overlap of the extended object between images, so it cannot be removed from the flat.
enough that the sky flat is too noisy to be useful.
Producing a sky flat is much like producing any other flat. The images must have bias and dark current subtracted (and overscan if it is being used) then combined, rescaling each image to take into account different levels of background illumination.
It is important to scale the median of each image to the same value instead of scaling the mean because the presence of bright sources will affect the mean much more than the median.
One down side of producing sky flats is the need to process the science images twice. The first time all of the usual calibration steps except flat fielding are done, then the flats are produced, then each science image is flat corrected.
The partially reduced images are saved in a different folder than the completely reduced science images that were processed earlier.
The images for this example were taken the same night as the other images in "Example 2" in earlier notebooks.
First, we set up some of the locations we will need.
ex2_calibrated = Path('example2-reduced')
sky_flat_bad_raw = Path('sky_flat_good_raw')
sky_flat_bad_working = Path('sky_flat_good_working')
sky_flat_bad_working.mkdir(exist_ok=True)
Next, load the combined bias and combined dark for this night. Recall that the combined dark for this night was bias-subtracted because it needed to be scaled for the flat images (see this notebook for more detail).
All of the science exposures this night had the same exposure time, 90 sec.
combined_bias = CCDData.read(ex2_calibrated / 'combined_bias.fit')
combined_dark = CCDData.read(ex2_calibrated / 'combined_dark_90.000.fit')
The telescope tracking changed during this night. Tracking was excellent for observations of Kelt 16b, making the images terrible for sky flats, but excellent for illustrating the failure of sky flats under some circumstances.
ifc_raw = ccdp.ImageFileCollection(sky_flat_bad_raw)
for ccd, name in ifc_raw.ccds(imagetyp='light', object='wasp 10 b', filter="r", return_fname=True):
reduced = ccdp.trim_image(ccd[:, :4096])
reduced = ccdp.subtract_bias(reduced, combined_bias)
reduced = ccdp.subtract_dark(reduced, combined_dark, exposure_time='exposure', exposure_unit=u.second)
thresh = detect_threshold(an_im, 2)
segm = detect_sources(an_im, thresh, 30)
reduced.data[segm.data > 0] = np.nan
reduced.data = reduced.data.astype('float32')
reduced.write(sky_flat_bad_working / name)
The combination settings here are important. Either combine by averaging and sigma clip or combine by median. Either should ensure that stars do not show up in your final flat as long as there is enough offset between the images. Images need to be scaled so that the median is the same for each image. Typically, a value of one is chosen as the common value.
ifc_working = ccdp.ImageFileCollection(sky_flat_bad_working)
to_combine = [ccd for ccd in ifc_working.ccds()]
def inv_median(array):
return 1 / np.nanmedian(array)
sky_flat = ccdp.combine(to_combine, scale=inv_median,
sigma_clip=True, sigma_clip_low_thresh=3, sigma_clip_high_thresh=3,
sigma_clip_func=np.nanmedian, sigma_clip_dev_func=mad_std,
mem_limit=2e9
)
from convenience_functions import show_image
an_im = CCDData.read(sky_flat_bad_working / 'wasp-10-b-S001-R001-C050-r.fit')
show_image(an_im, cmap='gray')
from photutils import detect_threshold, detect_sources, source_properties
foo = detect_threshold(an_im, 2)
arf = detect_sources(an_im, foo, 30)
show_image(arf.data > 0, cmap='gray', is_mask=True)
(arf.data > 0).sum()
moo = source_properties(an_im.data, arf)
moo.to_table()
show_image(sky_flat, cmap='gray')
sky_flat.write('supposed_to_be_good_but_has_streaks.fits')