LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
cameraGeometry.py
Go to the documentation of this file.
1 # This file is part of jointcal.
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 """Code to convert jointcal's output WCS models to distortion maps that can be
23 used by afw CameraGeom.
24 """
25 import numpy as np
26 
27 from lsst.afw import cameraGeom
28 import lsst.afw.geom
29 import astshim as ast
30 import lsst.log
31 from lsst.geom import SpherePoint, Point2D, radians
32 
33 _LOG = lsst.log.Log.getLogger(__name__)
34 
35 
37  """Convert a jointcal `~lsst.afw.geom.SkyWcs` into a distortion model and
38  detector positions (TODO) that can be used by `~lsst.afw.cameraGeom`.
39 
40  Because this code only operates on the WCS, it is independent of the
41  format of the persisted output (e.g. gen2 separate files vs. gen3 bundled
42  visits).
43 
44  Parameters
45  ----------
46  wcsList : `list` [`lsst.afw.geom.SkyWcs`]
47  The WCS to use to compute the distortion model from, preferably from
48  multiple visits on the same tract.
49  detectors : `list` [`int`]
50  Detector ids that correspond one-to-one with ``wcsList``.
51  camera : `lsst.afw.cameraGeom.Camera`
52  The camera these WCS were fit for.
53  n : `int`
54  Number of points to compute the pixel scale at, along the +y axis.
55  """
56  def __init__(self, wcsList, detectors, camera, n=100):
57  self.wcsListwcsList = wcsList
58  self.cameracamera = camera
59  self.detectorsdetectors = detectors
60  self.maxFocalRadiusmaxFocalRadius = self.cameracamera.computeMaxFocalPlaneRadius()
61  self.nn = n
62  # the computed radius and pixel scales
63  self.fieldAnglefieldAngle = None # degrees
64  self.radialScaleradialScale = None # arcsec
65  self.tangentialScaletangentialScale = None # arcsec
66  # the computed values for every input wcs
67  self.fieldAnglesfieldAngles = None
68  self.radialScalesradialScales = None
69  self.tangentialScalestangentialScales = None
70  self.fieldAngleStdfieldAngleStd = None
71  self.radialScaleStdradialScaleStd = None
72  self.tangentialScaleStdtangentialScaleStd = None
73 
74  self.loglog = _LOG.getChild("CameraModel")
75 
77  """Calculate the afw cameraGeom distortion model to be included in an
78  on-disk camera model.
79 
80 
81  PLACEHOLDER: This may be as simple as running `computePixelScale` and
82  then doing a numpy polynomial fit to it for the cameraGeom input.
83  However, we need to check details of how that distortion model is
84  stored in a Camera.
85  e.g.: np.polyfit(self.fieldAngle, self.radialScale, poly_degree))
86  """
87  raise NotImplementedError("not yet!")
88 
89  def computePixelScale(self):
90  """Compute the radial and tangential pixel scale by averaging over
91  multiple jointcal WCS models.
92 
93  Also computes the standard deviation and logs any WCS that are
94  significant outliers.
95  The calculations are stored in the ``fieldAngle[s]``,
96  ``radialScale[s]``, and ``tangentialScale[s]`` member variables.
97  """
98  self.fieldAnglesfieldAngles = []
99  self.radialScalesradialScales = []
100  self.tangentialScalestangentialScales = []
101  for id, wcs in zip(self.detectorsdetectors, self.wcsListwcsList):
102  fieldAngle, radial, tangential = self._computeDetectorPixelScale_computeDetectorPixelScale(id, wcs)
103  self.fieldAnglesfieldAngles.append(fieldAngle)
104  self.radialScalesradialScales.append(radial)
105  self.tangentialScalestangentialScales.append(tangential)
106  # TODO: For now, don't worry about small differences in the computed
107  # field angles vs. their respective radial/tangential scales, just
108  # assume all fieldAngle positions are "close enough" and warn if not.
109  self.fieldAnglefieldAngle = np.mean(self.fieldAnglesfieldAngles, axis=0)
110  self.fieldAngleStdfieldAngleStd = np.std(self.fieldAnglesfieldAngles, axis=0)
111  if self.fieldAngleStdfieldAngleStd.max() > 1e-4:
112  self.loglog.warning("Large stddev in computed field angles between visits (max: %s degree).",
113  self.fieldAngleStdfieldAngleStd.max())
114  # import os; print(os.getpid()); import ipdb; ipdb.set_trace();
115  self.radialScaleradialScale = np.mean(self.radialScalesradialScales, axis=0)
116  self.radialScaleStdradialScaleStd = np.std(self.radialScalesradialScales, axis=0)
117  if self.radialScaleStdradialScaleStd.max() > 1e-4:
118  self.loglog.warning("Large stddev in computed radial scales between visits (max: %s arcsec).",
119  self.radialScaleStdradialScaleStd.max())
120  self.tangentialScaletangentialScale = np.mean(self.tangentialScalestangentialScales, axis=0)
121  self.tangentialScaleStdtangentialScaleStd = np.std(self.tangentialScalestangentialScales, axis=0)
122  if self.tangentialScaleStdtangentialScaleStd.max() > 1e-4:
123  self.loglog.warning("Large stddev in computed tangential scales between visits (max: %s arcsec).",
124  self.tangentialScaleStdtangentialScaleStd.max())
125 
126  def computeCameraPixelScale(self, detector_id=30):
127  """Compute the radial and tangential pixel scales using the distortion
128  model supplied with the camera.
129 
130  This is designed to be directly comparable with the results of
131  `~CameraModel.computePixelScale`.
132 
133  Parameters
134  ----------
135  detector_id: `int`
136  Detector identifier for the detector_id to use for the calculation.
137 
138  Returns
139  -------
140  fieldAngle : `numpy.ndarray`
141  Field angles in degrees.
142  radialScale : `numpy.ndarray`
143  Radial direction pixel scales in arcseconds/pixel.
144  tangentialScale : `numpy.ndarray`
145  Tangential direction pixel scales in arcseconds/pixel.
146  """
147  # Make a trivial SkyWcs to get a field angle->sky transform from.
148  iwcToSkyWcs = lsst.afw.geom.makeSkyWcs(Point2D(0, 0), SpherePoint(0, 0, radians),
149  lsst.afw.geom.makeCdMatrix(1.0 * radians, 0 * radians, True))
150  iwcToSkyMap = iwcToSkyWcs.getFrameDict().getMapping("PIXELS", "SKY")
151  skyFrame = iwcToSkyWcs.getFrameDict().getFrame("SKY")
152 
153  # Extract the transforms that are defined just on the camera.
154  pixSys = self.cameracamera[detector_id].makeCameraSys(cameraGeom.PIXELS)
155  pixelsToFocal = self.cameracamera.getTransform(pixSys, cameraGeom.FOCAL_PLANE)
156  focalToField = self.cameracamera.getTransform(cameraGeom.FOCAL_PLANE, cameraGeom.FIELD_ANGLE)
157 
158  # Build a SkyWcs that combines each of the above components.
159  pixelFrame = ast.Frame(2, "Domain=PIXELS")
160  focalFrame = ast.Frame(2, "Domain=FOCAL")
161  iwcFrame = ast.Frame(2, "Domain=IWC")
162  frameDict = ast.FrameDict(pixelFrame)
163  frameDict.addFrame("PIXELS", pixelsToFocal.getMapping(), focalFrame)
164  frameDict.addFrame("FOCAL", focalToField.getMapping(), iwcFrame)
165  frameDict.addFrame("IWC", iwcToSkyMap, skyFrame)
166  wcs = lsst.afw.geom.SkyWcs(frameDict)
167 
168  return self._computeDetectorPixelScale_computeDetectorPixelScale(detector_id, wcs)
169 
170  def _computeDetectorPixelScale(self, detector_id, wcs):
171  """Compute pixel scale in radial and tangential directions as a
172  function of field angle.
173 
174  Parameters
175  ----------
176  detector_id: `int`
177  Detector identifier for the detector of this wcs.
178  wcs : `lsst.afw.geom.SkyWcs`
179  Full focal-plane model to compute pixel scale on.
180 
181  Returns
182  -------
183  fieldAngle : `numpy.ndarray`
184  Field angles in degrees.
185  radialScale : `numpy.ndarray`
186  Radial direction pixel scales in arcseconds/pixel.
187  tangentialScale : `numpy.ndarray`
188  Tangential direction pixel scales in arcseconds/pixel.
189 
190  Notes
191  -----
192  Pixel scales are calculated from finite differences only along the +y
193  focal plane direction.
194  """
195  focalToSky = wcs.getFrameDict().getMapping('FOCAL', 'SKY')
196  mmPerPixel = self.cameracamera[detector_id].getPixelSize()
197 
198  focalToPixels = wcs.getFrameDict().getMapping('FOCAL', 'PIXELS')
199  trans = wcs.getTransform() # Pixels to Sky as Point2d -> SpherePoint
200  boresight = trans.applyForward(Point2D(focalToPixels.applyForward([0, 0])))
201 
202  rs = np.linspace(0, self.maxFocalRadiusmaxFocalRadius, self.nn) # focal plane units
203  fieldAngle = np.zeros_like(rs)
204  radialScale = np.zeros_like(rs)
205  tangentialScale = np.zeros_like(rs)
206  for i, r in enumerate(rs):
207  # point on the sky at position r along the focal plane +y axis
208  sp1 = SpherePoint(*focalToSky.applyForward(Point2D([0, r])), radians)
209  # point on the sky one pixel further along the focal plane +y axis
210  sp2 = SpherePoint(*focalToSky.applyForward(Point2D([0, r + mmPerPixel.getY()])), radians)
211  # point on the sky one pixel off of the focal plane +y axis at r
212  sp3 = SpherePoint(*focalToSky.applyForward(Point2D([mmPerPixel.getX(), r])), radians)
213  fieldAngle[i] = boresight.separation(sp1).asDegrees()
214  radialScale[i] = sp1.separation(sp2).asArcseconds()
215  tangentialScale[i] = sp1.separation(sp3).asArcseconds()
216  return fieldAngle, radialScale, tangentialScale
int max
A FrameSet whose frames can be referenced by domain name.
Definition: FrameDict.h:67
Frame is used to represent a coordinate system.
Definition: Frame.h:157
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
Definition: SkyWcs.h:117
Point in an unspecified spherical coordinate system.
Definition: SpherePoint.h:57
def computeCameraPixelScale(self, detector_id=30)
def __init__(self, wcsList, detectors, camera, n=100)
def _computeDetectorPixelScale(self, detector_id, wcs)
static Log getLogger(Log const &logger)
Definition: Log.h:772
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Definition: functional.cc:33
std::shared_ptr< SkyWcs > makeSkyWcs(daf::base::PropertySet &metadata, bool strip=false)
Construct a SkyWcs from FITS keywords.
Definition: SkyWcs.cc:521
Eigen::Matrix2d makeCdMatrix(lsst::geom::Angle const &scale, lsst::geom::Angle const &orientation=0 *lsst::geom::degrees, bool flipX=false)
Make a WCS CD matrix.
Definition: SkyWcs.cc:133
Point< double, 2 > Point2D
Definition: Point.h:324
Definition: Log.h:717