LSSTApplications  19.0.0-10-g920eed2,19.0.0-11-g48a0200+2,19.0.0-18-gfc4e62b+13,19.0.0-2-g3b2f90d+2,19.0.0-2-gd671419+5,19.0.0-20-g5a5a17ab+11,19.0.0-21-g2644856+13,19.0.0-23-g84eeccb+1,19.0.0-24-g878c510+1,19.0.0-25-g6c8df7140,19.0.0-25-gb330496+1,19.0.0-3-g2b32d65+5,19.0.0-3-g8227491+12,19.0.0-3-g9c54d0d+12,19.0.0-3-gca68e65+8,19.0.0-3-gcfc5f51+5,19.0.0-3-ge110943+11,19.0.0-3-ge74d124,19.0.0-3-gfe04aa6+13,19.0.0-30-g9c3fd16+1,19.0.0-4-g06f5963+5,19.0.0-4-g3d16501+13,19.0.0-4-g4a9c019+5,19.0.0-4-g5a8b323,19.0.0-4-g66397f0+1,19.0.0-4-g8278b9b+1,19.0.0-4-g8557e14,19.0.0-4-g8964aba+13,19.0.0-4-ge404a01+12,19.0.0-5-g40f3a5a,19.0.0-5-g4db63b3,19.0.0-5-gfb03ce7+13,19.0.0-6-gbaebbfb+12,19.0.0-61-gec4c6e08+1,19.0.0-7-g039c0b5+11,19.0.0-7-gbea9075+4,19.0.0-7-gc567de5+13,19.0.0-71-g41c0270,19.0.0-9-g2f02add+1,19.0.0-9-g463f923+12,w.2020.22
LSSTDataManagementBasePackage
crosstalk.py
Go to the documentation of this file.
1 #
2 # This file is part of obs_decam.
3 #
4 # Developed for the LSST Data Management System.
5 # This product includes software developed by the LSST Project
6 # (http://www.lsst.org).
7 # Copyright 2018 LSST Corporation.
8 #
9 # This program is free software: you can redistribute it and/or modify
10 # it under the terms of the GNU General Public License as published by
11 # the Free Software Foundation, either version 3 of the License, or
12 # (at your option) any later version.
13 #
14 # This program is distributed in the hope that it will be useful,
15 # but WITHOUT ANY WARRANTY; without even the implied warranty of
16 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 # GNU General Public License for more details.
18 #
19 # You should have received a copy of the GNU General Public License
20 # along with this program. If not, see <http://www.gnu.org/licenses/>.
21 #
22 
23 """Apply crosstalk corrections to DECam images.
24 
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)
27 """
28 
29 
30 from collections import defaultdict
31 from astropy.io import fits
32 import datetime as dt
33 import os
34 import argparse
35 import re
36 import numpy as np
37 
38 import lsst.log
39 import lsst.pipe.base as pipeBase
40 from lsst.ip.isr import CrosstalkConfig, CrosstalkTask, IsrTask
41 from lsst.obs.decam import DecamMapper
42 from lsst.utils import getPackageDir
43 import lsst.afw.math as afwMath
44 
45 
46 __all__ = ["DecamCrosstalkTask"]
47 
48 
50  """Configuration for DECam crosstalk removal.
51  """
52 
53  def getSourcesAndCoeffsFile(self, filename='DECam_xtalk_20130606.txt'):
54  """File containing DECam crosstalk coefficients.
55 
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
59 
60  Parameters
61  ----------
62  filename : `str`, optional
63  File containing the decam crosstalk coefficients, from NOAO.
64 
65  Returns
66  -------
67  result : `str`
68  Full path to filename.
69  """
70  mapper = DecamMapper()
71  packageName = mapper.getPackageName()
72  packageDir = getPackageDir(packageName)
73  return os.path.join(packageDir, 'decam', filename)
74 
76  """Read crosstalk sources and coefficients from DECam-specific file.
77 
78  Returns
79  -------
80  sources : `defaultdict`
81  Sources of crosstalk read from crosstalk_file.
82  coeffs : `dict`
83  Linear crosstalk coefficients corresponding to sources.
84  """
85  crosstalk_file = self.getSourcesAndCoeffsFile()
86  sources = defaultdict(list)
87  coeffs = {}
88  log = lsst.log.Log.getLogger('obs.decam.DecamCrosstalkConfig')
89  log.info('Reading crosstalk coefficient data')
90  with open(crosstalk_file) as f:
91  for line in f:
92  li = line.strip()
93  if not li.startswith('#'):
94  elem = li.split()
95  # The xtalk file has image areas like 'ccd01A'; preserve only '01A'
96  sources[elem[0][3:]].append(elem[1][3:])
97  # The linear crosstalk coefficients
98  coeffs[(elem[0][3:], elem[1][3:])] = float(elem[2])
99  return sources, coeffs
100 
101 
103  """Remove crosstalk from a raw DECam image.
104 
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.
108  """
109  ConfigClass = DecamCrosstalkConfig
110  _DefaultName = 'decamCrosstalk'
111 
112  def prepCrosstalk(self, dataRef):
113  """Retrieve source image data and coefficients to prepare for crosstalk correction.
114 
115  This is called by readIsrData. The crosstalkSources land in isrData.
116 
117  Parameters
118  ----------
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.
122 
123  Returns
124  -------
125  crosstalkSources : `defaultdict`
126  Contains image data and corresponding crosstalk coefficients for
127  each crosstalk source of the given dataRef victim.
128  """
129  # Retrieve crosstalk sources and coefficients for the whole focal plane
130  sources, coeffs = self.config.getSourcesAndCoeffs()
131 
132  # Define a butler to prepare for accessing crosstalk 'source' image data
133  try:
134  butler = dataRef.getButler()
135  except RuntimeError:
136  self.log.fatal('Cannot get a Butler from the dataRef provided')
137 
138  # Retrieve visit and ccdnum from 'victim' so we can look up 'sources'
139  victim_exposure = dataRef.get()
140  det = victim_exposure.getDetector()
141  visit = dataRef.dataId['visit']
142  ccdnum = det.getId()
143  ccdnum_str = '%02d' % ccdnum
144  crosstalkSources = defaultdict(list)
145  decamisr = IsrTask() # needed for source_exposure overscan correction
146  for amp in det:
147  victim = ccdnum_str + amp.getName()
148  for source in sources[victim]:
149  source_ccdnum = int(source[:2])
150  try:
151  source_exposure = butler.get('raw', dataId={'visit': visit, 'ccdnum': source_ccdnum})
152  except RuntimeError:
153  self.log.warn('Cannot access source %d, SKIPPING IT' % source_ccdnum)
154  else:
155  self.log.info('Correcting victim %s from crosstalk source %s' % (victim, source))
156  if 'A' in source:
157  amp_idx = 0
158  elif 'B' in source:
159  amp_idx = 1
160  else:
161  self.log.fatal('DECam source amp name does not contain A or B, cannot proceed')
162  source_amp = source_exposure.getDetector()[amp_idx]
163  # If doCrosstalkBeforeAssemble=True, then use getRawBBox(). Otherwise, getRawDataBBox().
164  source_dataBBox = source_amp.getRawBBox()
165  # Check to see if source overscan has been corrected, and if not, correct it first
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'
169  % source)
170  source_image = source_exposure.getMaskedImage().getImage()
171 
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
175 
176  @pipeBase.timeMethod
177  def run(self, exposure, crosstalkSources=None):
178  """Perform crosstalk correction on a DECam exposure and its corresponding dataRef.
179 
180  Parameters
181  ----------
182  exposure : `lsst.afw.image.Exposure`
183  Exposure to correct.
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
190  by prepCrosstalk.
191 
192  Returns
193  -------
194  `lsst.pipe.base.Struct`
195  Struct with components:
196  - ``exposure``: The exposure after crosstalk correction has been
197  applied (`lsst.afw.image.Exposure`).
198  """
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."
203  # Load data from crosstalk 'victim' exposure we want to correct
204  det = exposure.getDetector()
205  for amp in det:
206  ccdnum = det.getId()
207  ccdnum_str = '%02d' % ccdnum
208  victim = ccdnum_str + amp.getName()
209  # If doCrosstalkBeforeAssemble=True, then use getRawBBox(). Otherwise, getBBox().
210  dataBBox = amp.getRawBBox()
211  # Check to see if victim overscan has been corrected, and if not, correct it first
212  if not exposure.getMetadata().exists('OVERSCAN'):
213  decamisr = IsrTask()
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)
219  # Load data from crosstalk 'source' exposures
220  for source in crosstalkSources[victim]:
221  source_data, source_coeff, source_idx = source
222  victim_idx = 0
223  if 'A' in victim:
224  victim_idx = 0
225  elif 'B' in victim:
226  victim_idx = 1
227  else:
228  self.log.fatal('DECam victim amp name does not contain A or B, cannot proceed')
229 
230  if source_idx != victim_idx:
231  # Occurs with amp A and B mismatch; need to flip horizontally
232  source_data = afwMath.flipImage(source_data, flipLR=True, flipTB=False)
233  # Perform the linear crosstalk correction
234  try:
235  source_data *= source_coeff
236  victim_data -= source_data
237  except RuntimeError:
238  self.log.fatal('Crosstalk correction failed for victim %s from source %s'
239  % (victim, source))
240  return pipeBase.Struct(
241  exposure=exposure,
242  )
243 
244 
245 class DecamCrosstalkIO(pipeBase.Task):
246  """Remove crosstalk from a raw DECam Multi-Extension FITS file.
247 
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.
251  """
252  ConfigClass = DecamCrosstalkConfig
253  _DefaultName = 'decamCrosstalkIO'
254 
256  """Parse command-line arguments to run the crosstalk correction alone.
257 
258  Examples
259  --------
260  runCrosstalkAlone.py -f decam0411673.fits.fz -o decam0411673_xtalk_corrected.fits.fz
261  """
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}
269  return parsed
270 
271  def runCrosstalkAlone(self):
272  """Utility for crosstalk correction directly on a DECam MEF file.
273  """
274  log = lsst.log.Log.getLogger('obs.decam.DecamCrosstalkIO')
275  parsed = self.parseCrosstalkIOArgs()
276  filename = parsed['filename']
277  outfile = parsed['outfile']
278  mef = fits.open(filename)
279  sources, coeffs = self.config.getSourcesAndCoeffs()
280  corrected = subtractCrosstalkIO(mef, sources, coeffs)
281 
282  if outfile:
283  corrected.writeto(outfile, overwrite=True)
284  else:
285  log.warning('No outfile specified, not writing corrected image to disk')
286 
287 
288 def subtractCrosstalkIO(mef, sources, coeffs):
289  """Remove crosstalk from a DECam Multi-Extension FITS raw image.
290 
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.
294 
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.**
298 
299  Parameters
300  ----------
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.
305  coeffs : `list`
306  Linear crosstalk coefficients.
307 
308  Returns
309  -------
310  mef : `astropy.io.fits.hdu.hdulist.HDUList`
311  New MEF image with crosstalk corrected data and updated header.
312  """
313  log = lsst.log.Log.getLogger('obs.decam.subtractCrosstalkIO')
314  lowx = {} # Lower x pixel index of given chip region
315  highx = {} # Upper x pixel index of given chip region
316  extension = {} # FITS extension of given chip position from CCDNUM header
317  corrected = [0] * 63 # placeholder for corrected data
318 
319  for i in range(1, 63):
320  ccd = '%02d' % mef[i].header['CCDNUM']
321  extension[ccd] = i
322  for section in 'AB':
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)]:
329  for section in 'AB':
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]:
335  # If regions were read out on different sides of the CCD
336  # (i.e., an A and B mismatch) flip data in the x-direction
337  source_data = np.fliplr(source_data)
338  # Perform linear crosstalk correction
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))
341 
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])
346 
347  return mef
lsst.obs.decam.crosstalk.DecamCrosstalkConfig
Definition: crosstalk.py:49
lsst::log.log.logContinued.warn
def warn(fmt, *args)
Definition: logContinued.py:202
lsst.obs.decam.crosstalk.subtractCrosstalkIO
def subtractCrosstalkIO(mef, sources, coeffs)
Definition: crosstalk.py:288
lsst::log.log.logContinued.info
def info(fmt, *args)
Definition: logContinued.py:198
lsst::ip::isr.crosstalk.CrosstalkTask
Definition: crosstalk.py:131
lsst.obs.decam.decamMapper.DecamMapper
Definition: decamMapper.py:44
lsst.obs.decam.crosstalk.DecamCrosstalkIO
Definition: crosstalk.py:245
pex.config.history.format
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
ast::append
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Definition: functional.cc:33
lsst::log::Log::getLogger
static Log getLogger(Log const &logger)
Definition: Log.h:760
lsst.obs.decam.crosstalk.DecamCrosstalkTask.prepCrosstalk
def prepCrosstalk(self, dataRef)
Definition: crosstalk.py:112
lsst::utils::getPackageDir
std::string getPackageDir(std::string const &packageName)
return the root directory of a setup package
Definition: packaging.cc:33
lsst::ip::isr.isrTask.IsrTask
Definition: isrTask.py:785
lsst.obs.decam.crosstalk.DecamCrosstalkTask.run
def run(self, exposure, crosstalkSources=None)
Definition: crosstalk.py:177
lsst::log.log.logContinued.fatal
def fatal(fmt, *args)
Definition: logContinued.py:214
lsst.pipe.base.task.Task.config
config
Definition: task.py:149
lsst.pipe.base.task.Task.log
log
Definition: task.py:148
lsst::log
Definition: Log.h:706
lsst::utils
Definition: Backtrace.h:29
lsst.obs.decam.crosstalk.DecamCrosstalkTask
Definition: crosstalk.py:102
lsst::ip::isr
Definition: applyLookupTable.h:34
lsst::afw::math
Definition: statistics.dox:6
lsst.obs.decam
Definition: __init__.py:1
lsst.obs.decam.crosstalk.DecamCrosstalkConfig.getSourcesAndCoeffs
def getSourcesAndCoeffs(self)
Definition: crosstalk.py:75
lsst.obs.decam.crosstalk.DecamCrosstalkIO.runCrosstalkAlone
def runCrosstalkAlone(self)
Definition: crosstalk.py:271
lsst.pipe.base
Definition: __init__.py:1
lsst.obs.decam.crosstalk.DecamCrosstalkIO.parseCrosstalkIOArgs
def parseCrosstalkIOArgs(self)
Definition: crosstalk.py:255
lsst::ip::isr.crosstalk.CrosstalkConfig
Definition: crosstalk.py:37
lsst.obs.decam.crosstalk.DecamCrosstalkConfig.getSourcesAndCoeffsFile
def getSourcesAndCoeffsFile(self, filename='DECam_xtalk_20130606.txt')
Definition: crosstalk.py:53
lsst::afw::math::flipImage
std::shared_ptr< ImageT > flipImage(ImageT const &inImage, bool flipLR, bool flipTB)
Flip an image left–right and/or top–bottom.
Definition: rotateImage.cc:92