LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
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
23used by afw CameraGeom.
24"""
25__all__ = ["CameraModel"]
26
27import logging
28import numpy as np
29
30from lsst.afw import cameraGeom
31import lsst.afw.geom
32import astshim as ast
33from lsst.geom import SpherePoint, Point2D, radians
34
35_LOG = logging.getLogger(__name__)
36
37
39 """Convert a jointcal `~lsst.afw.geom.SkyWcs` into a distortion model and
40 detector positions (TODO) that can be used by `~lsst.afw.cameraGeom`.
41
42 Because this code only operates on the WCS, it is independent of the
43 format of the persisted output (e.g. gen2 separate files vs. gen3 bundled
44 visits).
45
46 Parameters
47 ----------
48 wcsList : `list` [`lsst.afw.geom.SkyWcs`]
49 The WCS to use to compute the distortion model from, preferably from
50 multiple visits on the same tract.
51 detectors : `list` [`int`]
52 Detector ids that correspond one-to-one with ``wcsList``.
53 camera : `lsst.afw.cameraGeom.Camera`
54 The camera these WCS were fit for.
55 n : `int`
56 Number of points to compute the pixel scale at, along the +y axis.
57 """
58 def __init__(self, wcsList, detectors, camera, n=100):
59 self.wcsList = wcsList
60 self.camera = camera
61 self.detectors = detectors
62 self.maxFocalRadius = self.camera.computeMaxFocalPlaneRadius()
63 self.n = n
64 # the computed radius and pixel scales
65 self.fieldAngle = None # degrees
66 self.radialScale = None # arcsec
67 self.tangentialScale = None # arcsec
68 # the computed values for every input wcs
69 self.fieldAngles = None
70 self.radialScales = None
71 self.tangentialScales = None
72 self.fieldAngleStd = None
73 self.radialScaleStd = None
75
76 self.log = _LOG.getChild("CameraModel")
77
79 """Calculate the afw cameraGeom distortion model to be included in an
80 on-disk camera model.
81
82 PLACEHOLDER: This may be as simple as running `computePixelScale` and
83 then doing a numpy polynomial fit to it for the cameraGeom input.
84 However, we need to check details of how that distortion model is
85 stored in a Camera. e.g.:
86 np.polyfit(self.fieldAngle, self.radialScale, poly_degree)
87 """
88 raise NotImplementedError("not yet!")
89
91 """Compute the radial and tangential pixel scale by averaging over
92 multiple jointcal WCS models.
93
94 Also computes the standard deviation and logs any WCS that are
95 significant outliers.
96 The calculations are stored in the ``fieldAngle[s]``,
97 ``radialScale[s]``, and ``tangentialScale[s]`` member variables.
98 """
99 self.fieldAngles = []
100 self.radialScales = []
101 self.tangentialScales = []
102 for id, wcs in zip(self.detectors, self.wcsList):
103 fieldAngle, radial, tangential = self._computeDetectorPixelScale(id, wcs)
104 self.fieldAngles.append(fieldAngle)
105 self.radialScales.append(radial)
106 self.tangentialScales.append(tangential)
107 # TODO: For now, don't worry about small differences in the computed
108 # field angles vs. their respective radial/tangential scales, just
109 # assume all fieldAngle positions are "close enough" and warn if not.
110 self.fieldAngle = np.mean(self.fieldAngles, axis=0)
111 self.fieldAngleStd = np.std(self.fieldAngles, axis=0)
112 if self.fieldAngleStd.max() > 1e-4:
113 self.log.warning("Large stddev in computed field angles between visits (max: %s degree).",
114 self.fieldAngleStd.max())
115 # import os; print(os.getpid()); import ipdb; ipdb.set_trace();
116 self.radialScale = np.mean(self.radialScales, axis=0)
117 self.radialScaleStd = np.std(self.radialScales, axis=0)
118 if self.radialScaleStd.max() > 1e-4:
119 self.log.warning("Large stddev in computed radial scales between visits (max: %s arcsec).",
120 self.radialScaleStd.max())
121 self.tangentialScale = np.mean(self.tangentialScales, axis=0)
122 self.tangentialScaleStd = np.std(self.tangentialScales, axis=0)
123 if self.tangentialScaleStd.max() > 1e-4:
124 self.log.warning("Large stddev in computed tangential scales between visits (max: %s arcsec).",
125 self.tangentialScaleStd.max())
126
127 def computeCameraPixelScale(self, detector_id=30):
128 """Compute the radial and tangential pixel scales using the distortion
129 model supplied with the camera.
130
131 This is designed to be directly comparable with the results of
132 `~CameraModel.computePixelScale`.
133
134 Parameters
135 ----------
136 detector_id: `int`
137 Detector identifier for the detector_id to use for the calculation.
138
139 Returns
140 -------
141 fieldAngle : `numpy.ndarray`
142 Field angles in degrees.
143 radialScale : `numpy.ndarray`
144 Radial direction pixel scales in arcseconds/pixel.
145 tangentialScale : `numpy.ndarray`
146 Tangential direction pixel scales in arcseconds/pixel.
147 """
148 # Make a trivial SkyWcs to get a field angle->sky transform from.
149 iwcToSkyWcs = lsst.afw.geom.makeSkyWcs(Point2D(0, 0), SpherePoint(0, 0, radians),
150 lsst.afw.geom.makeCdMatrix(1.0 * radians, 0 * radians, True))
151 iwcToSkyMap = iwcToSkyWcs.getFrameDict().getMapping("PIXELS", "SKY")
152 skyFrame = iwcToSkyWcs.getFrameDict().getFrame("SKY")
153
154 # Extract the transforms that are defined just on the camera.
155 pixSys = self.camera[detector_id].makeCameraSys(cameraGeom.PIXELS)
156 pixelsToFocal = self.camera.getTransform(pixSys, cameraGeom.FOCAL_PLANE)
157 focalToField = self.camera.getTransform(cameraGeom.FOCAL_PLANE, cameraGeom.FIELD_ANGLE)
158
159 # Build a SkyWcs that combines each of the above components.
160 pixelFrame = ast.Frame(2, "Domain=PIXELS")
161 focalFrame = ast.Frame(2, "Domain=FOCAL")
162 iwcFrame = ast.Frame(2, "Domain=IWC")
163 frameDict = ast.FrameDict(pixelFrame)
164 frameDict.addFrame("PIXELS", pixelsToFocal.getMapping(), focalFrame)
165 frameDict.addFrame("FOCAL", focalToField.getMapping(), iwcFrame)
166 frameDict.addFrame("IWC", iwcToSkyMap, skyFrame)
167 wcs = lsst.afw.geom.SkyWcs(frameDict)
168
169 return self._computeDetectorPixelScale(detector_id, wcs)
170
171 def _computeDetectorPixelScale(self, detector_id, wcs):
172 """Compute pixel scale in radial and tangential directions as a
173 function of field angle.
174
175 Parameters
176 ----------
177 detector_id: `int`
178 Detector identifier for the detector of this wcs.
179 wcs : `lsst.afw.geom.SkyWcs`
180 Full focal-plane model to compute pixel scale on.
181
182 Returns
183 -------
184 fieldAngle : `numpy.ndarray`
185 Field angles in degrees.
186 radialScale : `numpy.ndarray`
187 Radial direction pixel scales in arcseconds/pixel.
188 tangentialScale : `numpy.ndarray`
189 Tangential direction pixel scales in arcseconds/pixel.
190
191 Notes
192 -----
193 Pixel scales are calculated from finite differences only along the +y
194 focal plane direction.
195 """
196 focalToSky = wcs.getFrameDict().getMapping('FOCAL', 'SKY')
197 mmPerPixel = self.camera[detector_id].getPixelSize()
198
199 focalToPixels = wcs.getFrameDict().getMapping('FOCAL', 'PIXELS')
200 trans = wcs.getTransform() # Pixels to Sky as Point2d -> SpherePoint
201 boresight = trans.applyForward(Point2D(focalToPixels.applyForward([0, 0])))
202
203 rs = np.linspace(0, self.maxFocalRadius, self.n) # focal plane units
204 fieldAngle = np.zeros_like(rs)
205 radialScale = np.zeros_like(rs)
206 tangentialScale = np.zeros_like(rs)
207 for i, r in enumerate(rs):
208 # point on the sky at position r along the focal plane +y axis
209 sp1 = SpherePoint(*focalToSky.applyForward(Point2D([0, r])), radians)
210 # point on the sky one pixel further along the focal plane +y axis
211 sp2 = SpherePoint(*focalToSky.applyForward(Point2D([0, r + mmPerPixel.getY()])), radians)
212 # point on the sky one pixel off of the focal plane +y axis at r
213 sp3 = SpherePoint(*focalToSky.applyForward(Point2D([mmPerPixel.getX(), r])), radians)
214 fieldAngle[i] = boresight.separation(sp1).asDegrees()
215 radialScale[i] = sp1.separation(sp2).asArcseconds()
216 tangentialScale[i] = sp1.separation(sp3).asArcseconds()
217 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
_computeDetectorPixelScale(self, detector_id, wcs)
__init__(self, wcsList, detectors, camera, n=100)
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