LSSTApplications  17.0+11,17.0+34,17.0+56,17.0+57,17.0+59,17.0+7,17.0-1-g377950a+33,17.0.1-1-g114240f+2,17.0.1-1-g4d4fbc4+28,17.0.1-1-g55520dc+49,17.0.1-1-g5f4ed7e+52,17.0.1-1-g6dd7d69+17,17.0.1-1-g8de6c91+11,17.0.1-1-gb9095d2+7,17.0.1-1-ge9fec5e+5,17.0.1-1-gf4e0155+55,17.0.1-1-gfc65f5f+50,17.0.1-1-gfc6fb1f+20,17.0.1-10-g87f9f3f+1,17.0.1-11-ge9de802+16,17.0.1-16-ga14f7d5c+4,17.0.1-17-gc79d625+1,17.0.1-17-gdae4c4a+8,17.0.1-2-g26618f5+29,17.0.1-2-g54f2ebc+9,17.0.1-2-gf403422+1,17.0.1-20-g2ca2f74+6,17.0.1-23-gf3eadeb7+1,17.0.1-3-g7e86b59+39,17.0.1-3-gb5ca14a,17.0.1-3-gd08d533+40,17.0.1-30-g596af8797,17.0.1-4-g59d126d+4,17.0.1-4-gc69c472+5,17.0.1-6-g5afd9b9+4,17.0.1-7-g35889ee+1,17.0.1-7-gc7c8782+18,17.0.1-9-gc4bbfb2+3,w.2019.22
LSSTDataManagementBasePackage
fitsUtils.py
Go to the documentation of this file.
1 # This file is part of afw.
2 #
3 # Developed for the LSST Data Management System.
4 # This product includes software developed by the LSST Project
5 # (https://www.lsst.org).
6 # See the COPYRIGHT file at the top-level directory of this distribution
7 # for details of code ownership.
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 <https://www.gnu.org/licenses/>.
21 
22 __all__ = ['getByKey', 'setByKey', 'HeaderMap', 'HeaderAmpMap',
23  'HeaderDetectorMap', 'DetectorBuilder']
24 
25 import re
26 import warnings
27 
28 import lsst.geom
29 import lsst.afw.image as afwImage
30 import lsst.afw.table as afwTable
31 from lsst.afw.fits import readMetadata
32 import lsst.afw.cameraGeom as afwCameraGeom
33 import lsst.afw.geom
34 
35 
36 def getByKey(metadata, key):
37  """Wrapper for getting a value from a metadata object by key.
38 
39  Parameters
40  ----------
41  metadata : `lsst.daf.base.PropertySet`
42  Metadata object to query for value.
43  key : `str`
44  Key to use for value lookup.
45 
46  Returns
47  -------
48  value : `object`
49  Value associated with key, None if key does not exist.
50  """
51  mdKeys = metadata.paramNames()
52  if key in mdKeys:
53  return metadata.getScalar(key)
54  else:
55  return None
56 
57 
58 def setByKey(metadata, key, value, clobber):
59  """Wrapper for setting a value in a metadata object. Deals with case
60  where the key already exists.
61 
62  Parameters
63  ----------
64  metadata : `lsst.daf.base.PropertySet`
65  Metadata object ot modify in place.
66  key : `str`
67  Key to associate with value.
68  value : `object`
69  Value to assign in the metadata object.
70  clobber : `bool`
71  Clobber the value if the key already exists?
72  """
73  mdKeys = metadata.paramNames()
74  if key not in mdKeys or (key in mdKeys and clobber):
75  metadata.set(key, value)
76 
77 
78 class HeaderMap(dict):
79  """ Class to hold mapping of header cards to attributes"""
80 
81  def addEntry(self, keyname, attribute_name, default=None, transform=lambda x: x):
82  """Adds an entry to the registry.
83 
84  Parameters
85  ----------
86  keyname : `str`
87  Key used to retrieve the header record.
88  attribute_name : `str`
89  Name of the attribute to store the value in.
90  default : `object`
91  Default value to store if the header card is not available.
92  transform : `callable`
93  Transform to apply to the header value before assigning it to the
94  attribute.
95  """
96  self.__setitem__(attribute_name, {'keyName': keyname,
97  'default': default,
98  'transform': transform})
99 
100  def setAttributes(self, obj, metadata, doRaise=True):
101  """Sets the attributes on the given object given a metadata object.
102 
103  Parameters
104  ----------
105  obj : `object`
106  Object on which to operate in place.
107  metadata : `lsst.daf.base.PropertySet`
108  Metadata object used for applying the mapping.
109  doRaise : `bool`
110  Raise exceptions on calling methods on the input object that
111  do not exist?
112  """
113  for key, attrDict in self.items():
114  try:
115  value = getByKey(metadata, attrDict['keyName'])
116  # if attrDict['keyName'] == "RDCRNR" and value == 0:
117  # import ipdb; ipdb.set_trace()
118  if value is not None:
119  self._applyVal(obj, value, key, attrDict['transform'])
120  else:
121  # Only apply transform if the metadata has a value for this key
122  # otherwise assume the default value is transformed.
123  value = attrDict['default']
124  if value is not None:
125  self._applyVal(obj, value, key, lambda x: x)
126  except Exception as e:
127  if doRaise:
128  raise
129  else:
130  warnings.warn('WARNING: Failed to set %s attribute with %s value: %s' %
131  (key, value, str(e)))
132 
133  def _applyVal(self, obj, value, attrName, transform):
134  raise NotImplementedError('Must be implemented in sub-class')
135 
136 
138  """ Class to hold mapping of header cards to AmpInfoTable attributes
139  The amp info is stored using setters, thus calling the attribute as a function.
140  """
141 
142  def _applyVal(self, obj, value, attrName, transform):
143  getattr(obj, attrName)(transform(value))
144 
145 
147  """ Class to hold mapping of header cards to Detector attributes
148  Detector information is stored as attributes on a Config object.
149  """
150 
151  def _applyVal(self, obj, value, attrName, transform):
152  obj.__setattr__(attrName, transform(value))
153 
154 
156  """
157  Parameters
158  ----------
159  detectorFileName : `str`
160  FITS file containing the detector description.
161  May use [] notation to specify an extension in an MEF.
162  ampFileNameList : `list` of `str`
163  List of FITS file names to use in building the amps.
164  May contain duplicate entries if the raw data are assembled.
165  inAmpCoords : `bool`
166  True if raw data are in amp coordinates, False if raw data
167  are assembled into pseudo detector pixel arrays.
168  plateScale : `float`
169  Nominal platescale (arcsec/mm).
170  radialCoeffs : `iterable` of `float`
171  Radial distortion coefficients for a radial polynomial in
172  normalized units.
173  clobberMetadata : `bool`
174  Clobber metadata from input files if overridden in
175  _sanitizeMetadata().
176  doRaise : `bool`
177  Raise exception if not all non-defaulted keywords are defined?
178  """
179  def __init__(self, detectorFileName, ampFileNameList, inAmpCoords=True, plateScale=1.,
180  radialCoeffs=(0., 1.), clobberMetadata=False, doRaise=True):
181  self.inAmpCoords = inAmpCoords
184  self.detectorMetadata = readMetadata(detectorFileName)
186  self.detectorMetadata, clobber=clobberMetadata)
187  self.ampMetadataList = []
188  self.detector = None
189  self.doRaise = doRaise
190  for fileName in ampFileNameList:
191  self.ampMetadataList.append(readMetadata(fileName))
193  self.ampMetadataList[-1], clobber=clobberMetadata)
194  self.plateScale = plateScale
195  self.focalPlaneToField = self._makeRadialTransform(radialCoeffs)
196 
197  def _sanitizeHeaderMetadata(self, metadata, clobber):
198  """This method is called for all metadata and gives an opportunity to
199  add/modify header information for use downstream.
200 
201  Override this method if more than the default is needed.
202 
203  Parameters
204  ----------
205  metadata : `lsst.daf.base.PropertySet`
206  Metadata to read/modify
207  clobber : `bool`
208  Clobber keys that exist with default keys?
209  """
210  self._defaultSanitization(metadata, clobber)
211 
212  def _defaultSanitization(self, metadata, clobber):
213  """Does the default sanitization of the header metadata.
214 
215  Parameters
216  ----------
217  metadata : `lsst.daf.base.PropertySet`
218  Header metadata to extend/modify.
219  clobber : `bool`
220  Override values in existing header cards?
221  """
222 
223  if self.inAmpCoords:
224  # Deal with DTM to get flipX and flipY for assembly and add as 'FLIPX', 'FLIPY'
225  # The DTM array is a transformation matrix. As I understand it, it transforms between
226  # electronic and assembled coordintates. As such, a negative value in the DTM1_1 spot
227  # corresponds to a flip of the x-axis and a negative value in the DTM2_2 spot
228  # corresponds to a flip of the y-axis.
229  dtm1 = getByKey(metadata, 'DTM1_1')
230  dtm2 = getByKey(metadata, 'DTM2_2')
231  if dtm1 is not None and dtm2 is not None:
232  setByKey(metadata, 'FLIPX', dtm1 < 0, clobber)
233  setByKey(metadata, 'FLIPY', dtm2 < 0, clobber)
234  setByKey(metadata, 'RDCRNR', int(
235  afwTable.ReadoutCorner.LL), clobber)
236  else:
237  setByKey(metadata, 'FLIPX', False, clobber)
238  setByKey(metadata, 'FLIPY', True, clobber)
239  # I don't know how to figure out the read corner if already
240  # assembled
241  setByKey(metadata, 'RDCRNR', None, clobber)
242 
243  # Deal with NAXIS1, NAXIS2 to make rawBBox as 'RAWBBOX'
244  xext = getByKey(metadata, 'NAXIS1')
245  yext = getByKey(metadata, 'NAXIS2')
246  if xext is not None and yext is not None:
247  setByKey(metadata, 'RAWBBOX', '[%i:%i,%i:%i]'%(
248  1, xext, 1, yext), clobber)
249  # Deal with DTV1, DTV2 to make 'XYOFF
250  dtv1 = getByKey(metadata, 'DTV1')
251  dtv2 = getByKey(metadata, 'DTV2')
252  if dtv1 is not None and dtv2 is not None:
253  setByKey(metadata, 'XYOFF', [dtv1, dtv2], clobber)
254  # map biassec[1] to HOSCAN
255  # map biassec[3] to VOSCAN
256  # map biassec[2] to PRESCAN
257  if metadata.isArray('BIASSEC'):
258  keylist = ['HOSCAN', 'PRESCAN', 'VOSCAN']
259  biassecs = getByKey(metadata, 'BIASSEC')
260  for i, biassec in enumerate(biassecs):
261  setByKey(metadata, keylist[i], biassec, clobber)
262  else:
263  biassec = getByKey(metadata, 'BIASSEC')
264  if biassec is not None:
265  setByKey(metadata, 'HOSCAN', biassec, clobber)
266 
267  def _makeDefaultAmpMap(self):
268  """Make the default map from header information to amplifier
269  information.
270 
271  Returns
272  -------
273  headerAmMap : `HeaderAmpMap`
274  The HeaderAmpMap object containing the mapping
275  """
276  hMap = HeaderAmpMap()
277  emptyBBox = lsst.geom.BoxI()
278  mapList = [('EXTNAME', 'setName'),
279  ('DETSEC', 'setBBox', None, self._makeBbox),
280  ('GAIN', 'setGain', 1.),
281  ('RDNOISE', 'setReadNoise', 0.),
282  ('SATURATE', 'setSaturation', 2 << 15),
283  ('RDCRNR', 'setReadoutCorner', int(
284  afwTable.ReadoutCorner.LL), afwTable.ReadoutCorner),
285  ('LINCOEFF', 'setLinearityCoeffs', [0., 1.]),
286  ('LINTYPE', 'setLinearityType', 'POLY'),
287  ('RAWBBOX', 'setRawBBox', None, self._makeBbox),
288  ('DATASEC', 'setRawDataBBox', None, self._makeBbox),
289  ('FLIPX', 'setRawFlipX', False),
290  ('FLIPY', 'setRawFlipY', False),
291  ('XYOFF', 'setRawXYOffset', lsst.geom.ExtentI(0, 0), self._makeExt),
292  ('HOSCAN', 'setRawHorizontalOverscanBBox',
293  emptyBBox, self._makeBbox),
294  ('VOSCAN', 'setRawVerticalOverscanBBox',
295  emptyBBox, self._makeBbox),
296  ('PRESCAN', 'setRawPrescanBBox', emptyBBox, self._makeBbox),
297  ]
298  for tup in mapList:
299  hMap.addEntry(*tup)
300  return hMap
301 
302  def _makeDefaultDetectorMap(self):
303  """Make the default map from header information to detector information.
304 
305  Returns
306  -------
307  headerDetectorMap : `HeaderDetectorMap`
308  The HeaderDetectorMap object containing the mapping.
309  """
310  hMap = HeaderDetectorMap()
311  mapList = [('CCDNAME', 'name', 'ccdName'),
312  ('DETSIZE', 'bbox_x0', 0, self._getBboxX0),
313  ('DETSIZE', 'bbox_y0', 0, self._getBboxY0),
314  ('DETSIZE', 'bbox_x1', 0, self._getBboxX1),
315  ('DETSIZE', 'bbox_y1', 0, self._getBboxY1),
316  ('DETID', 'id', 0),
317  # DetectorConfig.detectorType is of type `int`, not
318  # `DetectorType`
319  ('OBSTYPE', 'detectorType', int(
320  afwCameraGeom.DetectorType.SCIENCE)),
321  ('SERSTR', 'serial', 'none'),
322  ('XPOS', 'offset_x', 0.),
323  ('YPOS', 'offset_y', 0.),
324  ('XPIX', 'refpos_x', 0.),
325  ('YPIX', 'refpos_y', 0.),
326  ('YAWDEG', 'yawDeg', 0.),
327  ('PITCHDEG', 'pitchDeg', 0.),
328  ('ROLLDEG', 'rollDeg', 0.),
329  ('XPIXSIZE', 'pixelSize_x', 1.),
330  ('YPIXSIZE', 'pixelSize_y', 1.),
331  ('TRNSPOSE', 'transposeDetector', False),
332  ]
333  for tup in mapList:
334  hMap.addEntry(*tup)
335  return hMap
336 
337  def _makeExt(self, extArr):
338  """Helper function to make an extent from an array
339 
340  Parameters
341  ----------
342  extArr : `array` of `int`
343  Length 2 array to use in creating the Extent object.
344 
345  Returns
346  -------
347  extent : `lsst.geom.Extent2I`
348  Extent constructed from the input list.
349  """
350  return lsst.geom.ExtentI(*extArr)
351 
352  def _makeBbox(self, boxString):
353  """Helper funtion to make a bounding box from a string representing a
354  FITS style bounding box.
355 
356  Parameters
357  ----------
358  boxString : `str`
359  String describing the bounding box.
360 
361  Returns
362  -------
363  bbox : `lsst.geom.Box2I`
364  The bounding box.
365  """
366  # strip off brackets and split into parts
367  x1, x2, y1, y2 = [int(el) for el in re.split(
368  '[:,]', boxString.strip()[1:-1])]
369  box = lsst.geom.BoxI(lsst.geom.PointI(x1, y1), lsst.geom.PointI(x2, y2))
370  # account for the difference between FITS convention and LSST convention for
371  # index of LLC.
372  box.shift(lsst.geom.Extent2I(-1, -1))
373  return box
374 
375  def _getBboxX0(self, boxString):
376  return self._makeBbox(boxString).getMinX()
377 
378  def _getBboxX1(self, boxString):
379  return self._makeBbox(boxString).getMaxX()
380 
381  def _getBboxY0(self, boxString):
382  return self._makeBbox(boxString).getMinY()
383 
384  def _getBboxY1(self, boxString):
385  return self._makeBbox(boxString).getMaxY()
386 
387  def _makeRadialTransform(self, radialCoeffs):
388  """Helper function to get the radial transform given the radial
389  polynomial coefficients given in the constructor.
390 
391  Parameters
392  ----------
393  radialCoeffs : `iterable` of `float`
394  List of coefficients describing a polynomial radial distortion in
395  normalized units. The first value must be 0.
396 
397  Returns
398  -------
399  transform : `lsst.afw.geom.TransformPoint2ToPoint2`
400  Transform object describing the radial distortion
401  """
402  pScaleRad = lsst.geom.arcsecToRad(self.plateScale)
403  return lsst.afw.geom.makeRadialTransform([el/pScaleRad for el in radialCoeffs])
404 
405  def buildDetector(self):
406  """Take all the information and build a Detector object.
407  The Detector object is necessary for doing things like assembly.
408 
409  Returns
410  -------
411  detector : `lsst.afw.cameraGeom.Detector`
412  The detector.
413  """
414  if self.detector is not None:
415  return self.detector
416 
417  schema = afwTable.AmpInfoTable.makeMinimalSchema()
418  ampInfo = afwTable.AmpInfoCatalog(schema)
419  for ampMetadata in self.ampMetadataList:
420  record = ampInfo.addNew()
421  self.defaultAmpMap.setAttributes(record, ampMetadata, self.doRaise)
422  record.setHasRawInfo(True)
423 
424  detConfig = afwCameraGeom.DetectorConfig()
425  self.defaultDetectorMap.setAttributes(
426  detConfig, self.detectorMetadata, self.doRaise)
427  self.detector = afwCameraGeom.makeDetector(
428  detConfig, ampInfo, self.focalPlaneToField)
429  return self.detector
430 
431  def makeExposure(self, im, mask=None, variance=None):
432  """Method for constructing an exposure object from an image and the
433  information contained in this class to construct the Detector.
434 
435  Parameters
436  ----------
437  im : `lsst.afw.image.Image`
438  Image used to construct the exposure.
439  mask : `lsst.afw.image.MaskU`
440  Optional mask plane.
441  variance : `lsst.afw.image.Image`
442  Optional variance plance as an image of the same type as im.
443 
444  Returns
445  -------
446  exposure : `lsst.afw.image.Exposure`
447  Constructed exposure (specific type will match that of ``im``).
448  """
449  if mask is None:
450  mask = afwImage.Mask(im.getDimensions())
451  if variance is None:
452  variance = im
453  mi = afwImage.makeMaskedImage(im, mask, variance)
454  detector = self.buildDetector()
455 
456  exp = afwImage.makeExposure(mi)
457  exp.setDetector(detector)
458  return exp
def makeExposure(self, im, mask=None, variance=None)
Definition: fitsUtils.py:431
def _makeRadialTransform(self, radialCoeffs)
Definition: fitsUtils.py:387
def addEntry(self, keyname, attribute_name, default=None, transform=lambda x:x)
Definition: fitsUtils.py:81
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Definition: functional.cc:33
MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > * makeMaskedImage(typename std::shared_ptr< Image< ImagePixelT >> image, typename std::shared_ptr< Mask< MaskPixelT >> mask=Mask< MaskPixelT >(), typename std::shared_ptr< Image< VariancePixelT >> variance=Image< VariancePixelT >())
A function to return a MaskedImage of the correct type (cf.
Definition: MaskedImage.h:1280
def _defaultSanitization(self, metadata, clobber)
Definition: fitsUtils.py:212
std::shared_ptr< Exposure< ImagePixelT, MaskPixelT, VariancePixelT > > makeExposure(MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > &mimage, std::shared_ptr< geom::SkyWcs const > wcs=std::shared_ptr< geom::SkyWcs const >())
A function to return an Exposure of the correct type (cf.
Definition: Exposure.h:457
def _sanitizeHeaderMetadata(self, metadata, clobber)
Definition: fitsUtils.py:197
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:78
def getByKey(metadata, key)
Definition: fitsUtils.py:36
std::shared_ptr< TransformPoint2ToPoint2 > makeRadialTransform(std::vector< double > const &forwardCoeffs, std::vector< double > const &inverseCoeffs)
A purely radial polynomial distortion.
def _applyVal(self, obj, value, attrName, transform)
Definition: fitsUtils.py:133
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects...
def setAttributes(self, obj, metadata, doRaise=True)
Definition: fitsUtils.py:100
constexpr double arcsecToRad(double x) noexcept
Definition: Angle.h:55
An integer coordinate rectangle.
Definition: Box.h:54
def __init__(self, detectorFileName, ampFileNameList, inAmpCoords=True, plateScale=1., radialCoeffs=(0., 1.), clobberMetadata=False, doRaise=True)
Definition: fitsUtils.py:180
def setByKey(metadata, key, value, clobber)
Definition: fitsUtils.py:58