22 """Code to convert jointcal's output WCS models to distortion maps that can be
23 used by afw CameraGeom.
31 from lsst.geom import SpherePoint, Point2D, radians
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`.
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
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.
54 Number of points to compute the pixel scale at, along the +y axis.
56 def __init__(self, wcsList, detectors, camera, n=100):
74 self.
loglog = _LOG.getChild(
"CameraModel")
77 """Calculate the afw cameraGeom distortion model to be included in an
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
85 e.g.: np.polyfit(self.fieldAngle, self.radialScale, poly_degree))
87 raise NotImplementedError(
"not yet!")
90 """Compute the radial and tangential pixel scale by averaging over
91 multiple jointcal WCS models.
93 Also computes the standard deviation and logs any WCS that are
95 The calculations are stored in the ``fieldAngle[s]``,
96 ``radialScale[s]``, and ``tangentialScale[s]`` member variables.
112 self.
loglog.
warning(
"Large stddev in computed field angles between visits (max: %s degree).",
118 self.
loglog.
warning(
"Large stddev in computed radial scales between visits (max: %s arcsec).",
123 self.
loglog.
warning(
"Large stddev in computed tangential scales between visits (max: %s arcsec).",
127 """Compute the radial and tangential pixel scales using the distortion
128 model supplied with the camera.
130 This is designed to be directly comparable with the results of
131 `~CameraModel.computePixelScale`.
136 Detector identifier for the detector_id to use for the calculation.
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.
150 iwcToSkyMap = iwcToSkyWcs.getFrameDict().getMapping(
"PIXELS",
"SKY")
151 skyFrame = iwcToSkyWcs.getFrameDict().
getFrame(
"SKY")
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)
159 pixelFrame =
ast.Frame(2,
"Domain=PIXELS")
160 focalFrame =
ast.Frame(2,
"Domain=FOCAL")
163 frameDict.addFrame(
"PIXELS", pixelsToFocal.getMapping(), focalFrame)
164 frameDict.addFrame(
"FOCAL", focalToField.getMapping(), iwcFrame)
165 frameDict.addFrame(
"IWC", iwcToSkyMap, skyFrame)
170 def _computeDetectorPixelScale(self, detector_id, wcs):
171 """Compute pixel scale in radial and tangential directions as a
172 function of field angle.
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.
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.
192 Pixel scales are calculated from finite differences only along the +y
193 focal plane direction.
195 focalToSky = wcs.getFrameDict().getMapping(
'FOCAL',
'SKY')
196 mmPerPixel = self.
cameracamera[detector_id].getPixelSize()
198 focalToPixels = wcs.getFrameDict().getMapping(
'FOCAL',
'PIXELS')
199 trans = wcs.getTransform()
200 boresight = trans.applyForward(
Point2D(focalToPixels.applyForward([0, 0])))
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):
210 sp2 =
SpherePoint(*focalToSky.applyForward(
Point2D([0, r + mmPerPixel.getY()])), 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
A FrameSet whose frames can be referenced by domain name.
Frame is used to represent a coordinate system.
A 2-dimensional celestial WCS that transform pixels to ICRS RA/Dec, using the LSST standard for pixel...
Point in an unspecified spherical coordinate system.
def computeCameraPixelScale(self, detector_id=30)
def computePixelScale(self)
def computeDistortionModel(self)
def __init__(self, wcsList, detectors, camera, n=100)
def _computeDetectorPixelScale(self, detector_id, wcs)
static Log getLogger(Log const &logger)
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
std::shared_ptr< SkyWcs > makeSkyWcs(daf::base::PropertySet &metadata, bool strip=false)
Construct a SkyWcs from FITS keywords.
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.
Point< double, 2 > Point2D