23 """Apply crosstalk corrections to DECam images.
25 Some of this code is based on DECam_crosstalk.py, used for the THELI pipeline
26 and written by Thomas Erben (private communication, 2018)
30 from collections
import defaultdict
31 from astropy.io
import fits
40 from lsst.ip.isr import CrosstalkConfig, CrosstalkTask, IsrTask
46 __all__ = [
"DecamCrosstalkTask"]
50 """Configuration for DECam crosstalk removal.
54 """File containing DECam crosstalk coefficients.
56 This text file is provided by NOAO in a particular format with
57 information about the DECam crosstalk coefficients. It is available at
58 http://www.ctio.noao.edu/noao/content/DECam-Calibration-Files
62 filename : `str`, optional
63 File containing the decam crosstalk coefficients, from NOAO.
68 Full path to filename.
71 packageName = mapper.getPackageName()
73 return os.path.join(packageDir,
'decam', filename)
76 """Read crosstalk sources and coefficients from DECam-specific file.
80 sources : `defaultdict`
81 Sources of crosstalk read from crosstalk_file.
83 Linear crosstalk coefficients corresponding to sources.
86 sources = defaultdict(list)
89 log.info(
'Reading crosstalk coefficient data')
90 with open(crosstalk_file)
as f:
93 if not li.startswith(
'#'):
96 sources[elem[0][3:]].
append(elem[1][3:])
98 coeffs[(elem[0][3:], elem[1][3:])] = float(elem[2])
99 return sources, coeffs
103 """Remove crosstalk from a raw DECam image.
105 This must be done after the overscan is corrected but before it is trimmed.
106 This version of crosstalk removal assumes you want to do it as part of a
107 regular stack-ISR workflow.
109 ConfigClass = DecamCrosstalkConfig
110 _DefaultName =
'decamCrosstalk'
113 """Retrieve source image data and coefficients to prepare for crosstalk correction.
115 This is called by readIsrData. The crosstalkSources land in isrData.
119 dataRef : `lsst.daf.persistence.butlerSubset.ButlerDataRef`
120 DataRef of exposure to correct which must include a dataId with
121 at least one visit and a ccdnum corresponding to the detector id.
125 crosstalkSources : `defaultdict`
126 Contains image data and corresponding crosstalk coefficients for
127 each crosstalk source of the given dataRef victim.
130 sources, coeffs = self.
config.getSourcesAndCoeffs()
134 butler = dataRef.getButler()
136 self.
log.
fatal(
'Cannot get a Butler from the dataRef provided')
139 victim_exposure = dataRef.get()
140 det = victim_exposure.getDetector()
141 visit = dataRef.dataId[
'visit']
143 ccdnum_str =
'%02d' % ccdnum
144 crosstalkSources = defaultdict(list)
147 victim = ccdnum_str + amp.getName()
148 for source
in sources[victim]:
149 source_ccdnum = int(source[:2])
151 source_exposure = butler.get(
'raw', dataId={
'visit': visit,
'ccdnum': source_ccdnum})
153 self.
log.
warn(
'Cannot access source %d, SKIPPING IT' % source_ccdnum)
155 self.
log.
info(
'Correcting victim %s from crosstalk source %s' % (victim, source))
161 self.
log.
fatal(
'DECam source amp name does not contain A or B, cannot proceed')
162 source_amp = source_exposure.getDetector()[amp_idx]
164 source_dataBBox = source_amp.getRawBBox()
166 if not source_exposure.getMetadata().exists(
'OVERSCAN'):
167 decamisr.overscanCorrection(source_exposure, source_amp)
168 self.
log.
info(
'Correcting source %s overscan before using to correct crosstalk'
170 source_image = source_exposure.getMaskedImage().getImage()
172 source_data = source_image.Factory(source_image, source_dataBBox, deep=
True)
173 crosstalkSources[victim].
append((source_data, coeffs[(victim, source)], amp_idx))
174 return crosstalkSources
177 def run(self, exposure, crosstalkSources=None):
178 """Perform crosstalk correction on a DECam exposure and its corresponding dataRef.
182 exposure : `lsst.afw.image.Exposure`
184 dataRef : `lsst.daf.persistence.butlerSubset.ButlerDataRef`
185 DataRef of exposure to correct which must include a dataId
186 with at least one visit and ccdnum.
187 crosstalkSources : `defaultdict`
188 Must contain image data and corresponding crosstalk coefficients for
189 each crosstalk source of the given dataRef victim. This is returned
194 `lsst.pipe.base.Struct`
195 Struct with components:
196 - ``exposure``: The exposure after crosstalk correction has been
197 applied (`lsst.afw.image.Exposure`).
199 self.
log.
info(
'Applying crosstalk correction')
200 assert crosstalkSources
is not None,
"Sources are required for DECam crosstalk correction; " \
201 "you must run CrosstalkTask via IsrTask which will " \
202 "call prepCrosstalk to get the sources."
204 det = exposure.getDetector()
207 ccdnum_str =
'%02d' % ccdnum
208 victim = ccdnum_str + amp.getName()
210 dataBBox = amp.getRawBBox()
212 if not exposure.getMetadata().exists(
'OVERSCAN'):
214 decamisr.overscanCorrection(exposure, amp)
215 self.
log.
warn(
'Overscan correction did not happen prior to crosstalk correction')
216 self.
log.
info(
'Correcting victim %s overscan before crosstalk' % victim)
217 image = exposure.getMaskedImage().getImage()
218 victim_data = image.Factory(image, dataBBox)
220 for source
in crosstalkSources[victim]:
221 source_data, source_coeff, source_idx = source
228 self.
log.
fatal(
'DECam victim amp name does not contain A or B, cannot proceed')
230 if source_idx != victim_idx:
235 source_data *= source_coeff
236 victim_data -= source_data
238 self.
log.
fatal(
'Crosstalk correction failed for victim %s from source %s'
240 return pipeBase.Struct(
246 """Remove crosstalk from a raw DECam Multi-Extension FITS file.
248 This must be done after the overscan is corrected but before it is trimmed.
249 NOTE: This version of crosstalk removal assumes you want to do it as an
250 standalone file operation, not as part of a regular stack-ISR workflow.
252 ConfigClass = DecamCrosstalkConfig
253 _DefaultName =
'decamCrosstalkIO'
256 """Parse command-line arguments to run the crosstalk correction alone.
260 runCrosstalkAlone.py -f decam0411673.fits.fz -o decam0411673_xtalk_corrected.fits.fz
262 parser = argparse.ArgumentParser()
263 parser.add_argument(
'-f',
'--filename',
264 help=
'Location on disk of the DECam MEF file you wish to crosstalk correct')
265 parser.add_argument(
'-o',
'--outfile',
266 help=
'Location on disk where crosstalk-corrected DECam MEF file will be saved')
267 args = parser.parse_args()
268 parsed = {
'filename': args.filename,
'outfile': args.outfile}
272 """Utility for crosstalk correction directly on a DECam MEF file.
276 filename = parsed[
'filename']
277 outfile = parsed[
'outfile']
278 mef = fits.open(filename)
279 sources, coeffs = self.config.getSourcesAndCoeffs()
283 corrected.writeto(outfile, overwrite=
True)
285 log.warning(
'No outfile specified, not writing corrected image to disk')
289 """Remove crosstalk from a DECam Multi-Extension FITS raw image.
291 This version uses astropy I/O rather than the Butler due to the cross-amp
292 crosstalk and ease of accessing Multi-Extension FITS data. It does not
293 integrate with isr and does not use the Butler.
295 WARNING: this will run fine on a non-overscan-corrected raw image, but the
296 crosstalk coefficients were computed post-overscan-correction. Therefore,
297 **the input raw image should be overscan corrected before this is run.**
301 mef : `astropy.io.fits.hdu.hdulist.HDUList`
302 One DECam Multi-Extension FITS image.
303 sources : `defaultdict`
304 Crosstalk source areas affecting flux in a certain victim area.
306 Linear crosstalk coefficients.
310 mef : `astropy.io.fits.hdu.hdulist.HDUList`
311 New MEF image with crosstalk corrected data and updated header.
319 for i
in range(1, 63):
320 ccd =
'%02d' % mef[i].header[
'CCDNUM']
323 ccdsec_str = mef[i].header[
'DATASEC' + section]
324 ccdsec_list = re.findall(
r"[\w']+", ccdsec_str)
325 lowx[ccd + section] = int(ccdsec_list[0]) - 1
326 highx[ccd + section] = int(ccdsec_list[1])
327 corrected[i] = mef[i].data.astype(np.float32)
328 for ccd
in [
'%02d' % i
for i
in range(1, 63)]:
330 victim = ccd + section
331 victim_data = corrected[extension[ccd]][:, lowx[victim]:highx[victim]]
332 for source
in sources[victim]:
333 source_data = mef[extension[source[:2]]].data[:, lowx[source]:highx[source]]
334 if lowx[victim] != lowx[source]:
337 source_data = np.fliplr(source_data)
339 victim_data[:, :] = victim_data - coeffs[(victim, source)] * source_data
340 log.info(
'Correcting victim %s from crosstalk source %s for HDU %s' % (victim, source, ccd))
342 for i
in range(1, 63):
343 mef[i].header[
'HISTORY'] =
'Crosstalk corrected on {0}'.
format(
344 dt.datetime.now().strftime(
"%Y-%m-%dT%H:%M:%S"))
345 mef[i] = fits.ImageHDU(header=mef[i].header, data=corrected[i])