LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
fitsUtils.py
Go to the documentation of this file.
1 import re, warnings
2 import lsst.afw.image as afwImage
3 import lsst.afw.table as afwTable
4 import lsst.afw.cameraGeom as afwCameraGeom
5 import lsst.afw.geom as afwGeom
6 
7 def getByKey(metadata, key):
8  """Wrapper for getting a value from a metadata object by key.
9  @param[in] metadata metadata object to query for value
10  @param[in] key key to use for value lookup
11  @returns value associated with key, None if key does not exist
12  """
13  mdKeys = metadata.paramNames()
14  if key in mdKeys:
15  return metadata.get(key)
16  else:
17  return None
18 
19 def setByKey(metadata, key, value, clobber):
20  """Wrapper for setting a value in a metadata object. Deals with case
21  where the key already exists.
22  @param[in, out] metadata metadata object ot modify in place.
23  @param[in] key key to associate with value
24  @param[in] value value to assign in the metadata object
25  @param[in] clobber Clobber the value if the key already exisists?
26  """
27  mdKeys = metadata.paramNames()
28  if key not in mdKeys or (key in mdKeys and clobber):
29  metadata.set(key, value)
30 
31 class HeaderMap(dict):
32  """ Class to hold mapping of header cards to attributes"""
33  def addEntry(self, keyname, attribute_name, default=None, transform=lambda x: x):
34  """Adds an entry to the registr
35  @param[in] keyname Key used to retrieve the header record
36  @param[in] attribute_name Name of the attribute to store the value in
37  @param[jn] default Default velue to store if the header card is not available
38  @param[in] transform Transform to apply to the header value before assigning it to the
39  attribute.
40  """
41  self.__setitem__(attribute_name, {'keyName':keyname,
42  'default':default,
43  'transform':transform})
44 
45  def setAttributes(self, obj, metadata, doRaise=True):
46  """Sets the attributes on the give object given a metadata object.
47  @param[in, out] obj Object on which to operate in place
48  @param[in] metadata Metadata object used for applying the mapping
49  @param[in] doRaise Raise exceptions on calling methods on the input object that do not exist?
50  """
51  for key, attrDict in self.iteritems():
52  try:
53  value = getByKey(metadata, attrDict['keyName'])
54  if value is not None:
55  self._applyVal(obj, value, key, attrDict['transform'])
56  else:
57  #Only apply transform if the metadata has a value for this key
58  #otherwise assume the default value is transformed.
59  value = attrDict['default']
60  self._applyVal(obj, value, key, lambda x: x)
61  except Exception, e:
62  if doRaise:
63  raise
64  else:
65  warnings.warn('WARNING: Failed to set %s attribute with %s value: %s'%
66  (key, value, e.message))
67 
68  def _applyVal(self, obj, value, attrName, transform):
69  raise NotImplementedError('Must be implemented in sub-class')
70 
72  """ Class to hold mapping of header cards to AmpInfoTable attributes
73  The amp info is stored using setters, thus calling the attribute as a function.
74  """
75  def _applyVal(self, obj, value, attrName, transform):
76  getattr(obj, attrName)(transform(value))
77 
79  """ Class to hold mapping of header cards to Detector attributes
80  Detector information is stored as attributes on a Config object.
81  """
82  def _applyVal(self, obj, value, attrName, transform):
83  obj.__setattr__(attrName, transform(value))
84 
85 class DetectorBuilder(object):
86  def __init__(self, detectorFileName, ampFileNameList, inAmpCoords=True, plateScale=1.,
87  radialCoeffs=(0., 1.), clobberMetadata=False, doRaise=True):
88  ''' @param[in] detectorFileName FITS file containing the detector description.
89  May use [] notation to specify an extension in an MEF.
90  @param[in] ampFileNameList List of FITS file names to use in building the amps.
91  May contain duplicate entries if the raw data are assembled.
92  @param[in] inAmpCoords Boolean, True if raw data are in amp coordinates, False if raw data
93  are assembled into pseudo detector pixel arrays
94  @param[in] plateScale Nominal platescale (arcsec/mm)
95  @param[in] radialCoeffs Radial distortion coefficients for a radial polynomial in normalized
96  units.
97  @param[in] clobberMetadata Clobber metadata from input files if overridden in the _sanitizeMetadata method
98  @param[in] doRaise Raise exception if not all non-defaulted keywords are defined? Default is True.
99  '''
100  self.inAmpCoords = inAmpCoords
103  self.detectorMetadata = afwImage.readMetadata(detectorFileName)
104  self._sanitizeHeaderMetadata(self.detectorMetadata, clobber=clobberMetadata)
105  self.ampMetadataList = []
106  self.detector = None
107  self.doRaise = doRaise
108  for fileName in ampFileNameList:
109  self.ampMetadataList.append(afwImage.readMetadata(fileName))
110  self._sanitizeHeaderMetadata(self.ampMetadataList[-1], clobber=clobberMetadata)
111  self.plateScale = plateScale
112  self.focalPlaneToPupil = self._makeRadialTransform(radialCoeffs)
113 
114  def _sanitizeHeaderMetadata(self, metadata, clobber):
115  """This method is called for all metadata and gives an opportunity to add/modify
116  header information for use downstream.
117  Override this method if more than the default is needed.
118  @param[in, out] metadata Metadata to read/modify
119  @param[in] clobber Clobber keys that exist with default keys?
120  """
121  self._defaultSanitization(metadata, clobber)
122 
123  def _defaultSanitization(self, metadata, clobber):
124  """Does the default sanitization of the header metadata.
125  @param[in,out] metadata Header metadata to extend/modify
126  @param[in] clobber Override values in existing header cards?
127  """
128 
129  if self.inAmpCoords:
130  #Deal with DTM to get flipX and flipY for assembly and add as 'FLIPX', 'FLIPY'
131  #The DTM array is a transformation matrix. As I understand it, it transforms between
132  #electronic and assembled coordintates. As such, a negative value in the DTM1_1 spot
133  #corresponds to a flip of the x-axis and a negative value in the DTM2_2 spot
134  #corresponds to a flip of the y-axis.
135  dtm1 = getByKey(metadata, 'DTM1_1')
136  dtm2 = getByKey(metadata, 'DTM2_2')
137  if dtm1 is not None and dtm2 is not None:
138  setByKey(metadata, 'FLIPX', dtm1 < 0, clobber)
139  setByKey(metadata, 'FLIPY', dtm2 < 0, clobber)
140  setByKey(metadata, 'RDCRNR', afwTable.LL, clobber)
141  else:
142  setByKey(metadata, 'FLIPX', False, clobber)
143  setByKey(metadata, 'FLIPY', True, clobber)
144  #I don't know how to figure out the read corner if already assembled
145  setByKey(metadata, 'RDCRNR', None, clobber)
146 
147  #Deal with NAXIS1, NAXIS2 to make rawBBox as 'RAWBBOX'
148  xext = getByKey(metadata, 'NAXIS1')
149  yext = getByKey(metadata, 'NAXIS2')
150  if xext is not None and yext is not None:
151  setByKey(metadata, 'RAWBBOX', '[%i:%i,%i:%i]'%(1, xext, 1, yext), clobber)
152  #Deal with DTV1, DTV2 to make 'XYOFF
153  dtv1 = getByKey(metadata, 'DTV1')
154  dtv2 = getByKey(metadata, 'DTV2')
155  if dtv1 is not None and dtv2 is not None:
156  setByKey(metadata, 'XYOFF', [dtv1, dtv2], clobber)
157  #map biassec[1] to HOSCAN
158  #map biassec[3] to VOSCAN
159  #map biassec[2] to PRESCAN
160  if metadata.isArray('BIASSEC'):
161  keylist = ['HOSCAN', 'PRESCAN', 'VOSCAN']
162  biassecs = getByKey(metadata, 'BIASSEC')
163  for i, biassec in enumerate(biassecs):
164  setByKey(metadata, keylist[i], biassec, clobber)
165  else:
166  biassec = getByKey(metadata, 'BIASSEC')
167  if biassec is not None:
168  setByKey(metadata, 'HOSCAN', biassec, clobber)
169 
171  """Make the default map from header information to amplifier information
172  @return The HeaderAmpMap object containing the mapping
173  """
174  hMap = HeaderAmpMap()
175  emptyBBox = afwGeom.BoxI()
176  mapList = [('EXTNAME', 'setName'),
177  ('DETSEC', 'setBBox', None, self._makeBbox),
178  ('GAIN', 'setGain', 1.),
179  ('RDNOISE', 'setReadNoise', 0.),
180  ('SATURATE', 'setSaturation', 2<<15),
181  ('RDCRNR', 'setReadoutCorner', afwTable.LL),
182  ('LINCOEFF', 'setLinearityCoeffs', [0., 1.]),
183  ('LINTYPE', 'setLinearityType', 'POLY'),
184  ('RAWBBOX', 'setRawBBox', None, self._makeBbox),
185  ('DATASEC', 'setRawDataBBox', None, self._makeBbox),
186  ('FLIPX', 'setRawFlipX', False),
187  ('FLIPY', 'setRawFlipY', False),
188  ('XYOFF', 'setRawXYOffset', afwGeom.ExtentI(0,0), self._makeExt),
189  ('HOSCAN', 'setRawHorizontalOverscanBBox', emptyBBox, self._makeBbox),
190  ('VOSCAN', 'setRawVerticalOverscanBBox', emptyBBox, self._makeBbox),
191  ('PRESCAN', 'setRawPrescanBBox', emptyBBox, self._makeBbox),
192  ]
193  for tup in mapList:
194  hMap.addEntry(*tup)
195  return hMap
196 
198  """Make the default map from header information to detector information
199  @return The HeaderDetectorMap object containing the mapping
200  """
201  hMap = HeaderDetectorMap()
202  mapList = [('CCDNAME', 'name', 'ccdName'),
203  ('DETSIZE', 'bbox_x0', 0, self._getBboxX0),
204  ('DETSIZE', 'bbox_y0', 0, self._getBboxY0),
205  ('DETSIZE', 'bbox_x1', 0, self._getBboxX1),
206  ('DETSIZE', 'bbox_y1', 0, self._getBboxY1),
207  ('DETID', 'id', 0),
208  ('OBSTYPE', 'detectorType', afwCameraGeom.SCIENCE),
209  ('SERSTR', 'serial', 'none'),
210  ('XPOS', 'offset_x', 0.),
211  ('YPOS', 'offset_y', 0.),
212  ('XPIX', 'refpos_x', 0.),
213  ('YPIX', 'refpos_y', 0.),
214  ('YAWDEG', 'yawDeg', 0.),
215  ('PITCHDEG', 'pitchDeg', 0.),
216  ('ROLLDEG', 'rollDeg', 0.),
217  ('XPIXSIZE', 'pixelSize_x', 1.),
218  ('YPIXSIZE', 'pixelSize_y', 1.),
219  ('TRNSPOSE', 'transposeDetector', False),
220  ]
221  for tup in mapList:
222  hMap.addEntry(*tup)
223  return hMap
224 
225  def _makeExt(self, extArr):
226  """Helper function to make an extent from an array
227  @param[in] extArr Length 2 array to use in creating the Extent object
228  @return Extent2I constructed from the input list
229  """
230  return afwGeom.ExtentI(*extArr)
231 
232  def _makeBbox(self, boxString):
233  """Helper funtion to make a bounding box from a string representing a FITS style bounding box
234  @param[in] boxString String describing the bounding box
235  @return Box2I for the bounding box
236  """
237  #strip off brackets and split into parts
238  x1, x2, y1, y2 = [int(el) for el in re.split('[:,]', boxString.strip()[1:-1])]
239  box = afwGeom.BoxI(afwGeom.PointI(x1, y1), afwGeom.PointI(x2, y2))
240  #account for the difference between FITS convention and LSST convention for
241  #index of LLC.
242  box.shift(afwGeom.Extent2I(-1, -1))
243  return box
244 
245 
246  def _getBboxX0(self, boxString):
247  return self._makeBbox(boxString).getMinX()
248 
249  def _getBboxX1(self, boxString):
250  return self._makeBbox(boxString).getMaxX()
251 
252  def _getBboxY0(self, boxString):
253  return self._makeBbox(boxString).getMinY()
254 
255  def _getBboxY1(self, boxString):
256  return self._makeBbox(boxString).getMaxY()
257 
258  def _makeRadialTransform(self, radialCoeffs):
259  """Helper function to get the radial transform given the radial polynomial coefficients given in
260  the constructor.
261  @param[in] radialCoeffs List of coefficients describing a polynomial radial distortion in
262  normalized units.
263  @return RadialXYTransform object describing the radial distortion
264  """
265  pScaleRad = afwGeom.arcsecToRad(self.plateScale)
266  return afwGeom.RadialXYTransform([el/pScaleRad for el in radialCoeffs])
267 
268  def buildDetector(self):
269  """Take all the information and build a Detector object. The Detector object is necessary for doing
270  things like assembly.
271  @return Detector object
272  """
273  if self.detector is not None:
274  return self.detector
275 
276  schema = afwTable.AmpInfoTable.makeMinimalSchema()
277  ampInfo = afwTable.AmpInfoCatalog(schema)
278  for ampMetadata in self.ampMetadataList:
279  record = ampInfo.addNew()
280  self.defaultAmpMap.setAttributes(record, ampMetadata, self.doRaise)
281  record.setHasRawInfo(True)
282 
283  detConfig = afwCameraGeom.DetectorConfig()
284  self.defaultDetectorMap.setAttributes(detConfig, self.detectorMetadata, self.doRaise)
285  self.detector = afwCameraGeom.makeDetector(detConfig, ampInfo, self.focalPlaneToPupil,
286  self.plateScale)
287  return self.detector
288 
289  def makeCalib(self):
290  """PLaceholder for subclasses to implement construction of a calib to associate with the exposure.
291  @return empty afwImage.Calib object
292  """
293  return afwImage.Calib()
294 
295  def makeExposure(self, im, mask=None, variance=None):
296  """Method for constructing an exposure object from an image and the information contained in this
297  class to construct the Detector and Calib objects.
298  @param[in] im Image used to construct the exposure
299  @param[in] mask Optional mask plane as a <askU
300  @param[in] variance Optional variance plance as an image of the same type as im
301  @param[out] Exposure object
302  """
303  if mask is None:
304  mask = afwImage.MaskU(im.getDimensions())
305  if variance is None:
306  variance = im
307  mi = afwImage.makeMaskedImage(im, mask, variance)
308  detector = self.buildDetector()
309 
311  calib = self.makeCalib()
312  exp = afwImage.makeExposure(mi, wcs)
313  exp.setCalib(calib)
314  exp.setDetector(detector)
315  return exp
MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > * makeMaskedImage(typename Image< ImagePixelT >::Ptr image, typename Mask< MaskPixelT >::Ptr mask=typename Mask< MaskPixelT >::Ptr(), typename Image< VariancePixelT >::Ptr variance=typename Image< VariancePixelT >::Ptr())
Definition: MaskedImage.h:1067
A custom container class for records, based on std::vector.
Definition: Catalog.h:94
An integer coordinate rectangle.
Definition: Box.h:53
Wcs::Ptr makeWcs(coord::Coord const &crval, geom::Point2D const &crpix, double CD11, double CD12, double CD21, double CD22)
Create a Wcs object from crval, crpix, CD, using CD elements (useful from python) ...
Definition: makeWcs.cc:141
boost::shared_ptr< daf::base::PropertySet > readMetadata(std::string const &fileName, int hdu=0, bool strip=false)
Return the metadata (header entries) from a FITS file.
A purely radial polynomial distortion, up to 6th order.
Definition: XYTransform.h:186
Exposure< ImagePixelT, MaskPixelT, VariancePixelT >::Ptr makeExposure(MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > &mimage, boost::shared_ptr< Wcs const > wcs=boost::shared_ptr< Wcs const >())
Definition: Exposure.h:308
double arcsecToRad(double x)
Definition: Angle.h:41