22 __all__ = [
"FitTanSipWcsTask",
"FitTanSipWcsConfig"]
33 from lsst.utils.timer
import timeMethod
34 from .setMatchDistance
import setMatchDistance
35 from .sip
import makeCreateWcsWithSip
39 """Config for FitTanSipWcsTask."""
40 order = pexConfig.RangeField(
41 doc=
"order of SIP polynomial",
46 numIter = pexConfig.RangeField(
47 doc=
"number of iterations of fitter (which fits X and Y separately, and so benefits from "
53 numRejIter = pexConfig.RangeField(
54 doc=
"number of rejection iterations",
59 rejSigma = pexConfig.RangeField(
60 doc=
"Number of standard deviations for clipping level",
65 maxScatterArcsec = pexConfig.RangeField(
66 doc=
"maximum median scatter of a WCS fit beyond which the fit fails (arcsec); "
67 "be generous, as this is only intended to catch catastrophic failures",
75 """Fit a TAN-SIP WCS given a list of reference object/source matches.
77 ConfigClass = FitTanSipWcsConfig
78 _DefaultName =
"fitWcs"
81 def fitWcs(self, matches, initWcs, bbox=None, refCat=None, sourceCat=None, exposure=None):
82 """Fit a TAN-SIP WCS from a list of reference object/source matches
86 matches : `list` of `lsst.afw.table.ReferenceMatch`
87 The following fields are read:
89 - match.first (reference object) coord
90 - match.second (source) centroid
92 The following fields are written:
94 - match.first (reference object) centroid,
95 - match.second (source) centroid
96 - match.distance (on sky separation, in radians)
98 initWcs : `lsst.afw.geom.SkyWcs`
100 bbox : `lsst.geom.Box2I`
101 the region over which the WCS will be valid (an lsst:afw::geom::Box2I);
102 if None or an empty box then computed from matches
103 refCat : `lsst.afw.table.SimpleCatalog`
104 reference object catalog, or None.
105 If provided then all centroids are updated with the new WCS,
106 otherwise only the centroids for ref objects in matches are updated.
107 Required fields are "centroid_x", "centroid_y", "coord_ra", and "coord_dec".
108 sourceCat : `lsst.afw.table.SourceCatalog`
109 source catalog, or None.
110 If provided then coords are updated with the new WCS;
111 otherwise only the coords for sources in matches are updated.
112 Required fields are "slot_Centroid_x", "slot_Centroid_y", and "coord_ra", and "coord_dec".
113 exposure : `lsst.afw.image.Exposure`
114 Ignored; present for consistency with FitSipDistortionTask.
118 result : `lsst.pipe.base.Struct`
119 with the following fields:
121 - ``wcs`` : the fit WCS (`lsst.afw.geom.SkyWcs`)
122 - ``scatterOnSky`` : median on-sky separation between reference
123 objects and sources in "matches" (`lsst.afw.geom.Angle`)
131 wcs = self.
initialWcsinitialWcs(matches, initWcs)
132 rejected = np.zeros(len(matches), dtype=bool)
133 for rej
in range(self.config.numRejIter):
134 sipObject = self.
_fitWcs_fitWcs([mm
for i, mm
in enumerate(matches)
if not rejected[i]], wcs)
135 wcs = sipObject.getNewWcs()
136 rejected = self.
rejectMatchesrejectMatches(matches, wcs, rejected)
137 if rejected.sum() == len(rejected):
138 raise RuntimeError(
"All matches rejected in iteration %d" % (rej + 1,))
140 "Iteration %d of astrometry fitting: rejected %f outliers, out of %d total matches.",
141 rej, rejected.sum(), len(rejected)
144 print(
"Plotting fit after rejection iteration %d/%d" % (rej + 1, self.config.numRejIter))
145 self.
plotFitplotFit(matches, wcs, rejected)
147 sipObject = self.
_fitWcs_fitWcs([mm
for i, mm
in enumerate(matches)
if not rejected[i]], wcs)
148 wcs = sipObject.getNewWcs()
150 print(
"Plotting final fit")
151 self.
plotFitplotFit(matches, wcs, rejected)
153 if refCat
is not None:
154 self.log.
debug(
"Updating centroids in refCat")
157 self.log.
warning(
"Updating reference object centroids in match list; refCat is None")
160 if sourceCat
is not None:
161 self.log.
debug(
"Updating coords in sourceCat")
164 self.log.
warning(
"Updating source coords in match list; sourceCat is None")
167 self.log.
debug(
"Updating distance in match list")
170 scatterOnSky = sipObject.getScatterOnSky()
172 if scatterOnSky.asArcseconds() > self.config.maxScatterArcsec:
173 raise pipeBase.TaskError(
174 "Fit failed: median scatter on sky = %0.3f arcsec > %0.3f config.maxScatterArcsec" %
175 (scatterOnSky.asArcseconds(), self.config.maxScatterArcsec))
177 return pipeBase.Struct(
179 scatterOnSky=scatterOnSky,
183 """Generate a guess Wcs from the astrometric matches
185 We create a Wcs anchored at the center of the matches, with the scale
186 of the input Wcs. This is necessary because matching returns only
187 matches with no estimated Wcs, and the input Wcs is a wild guess.
188 We're using the best of each: positions from the matches, and scale
193 matches : `list` of `lsst.afw.table.ReferenceMatch`
194 List of sources matched to references.
195 wcs : `lsst.afw.geom.SkyWcs`
200 newWcs : `lsst.afw.geom.SkyWcs`
201 Initial WCS guess from estimated crpix and crval.
207 crval += mm.first.getCoord().getVector()
208 crpix /= len(matches)
209 crval /= len(matches)
212 cdMatrix=wcs.getCdMatrix())
215 def _fitWcs(self, matches, wcs):
216 """Fit a Wcs based on the matches and a guess Wcs.
220 matches : `list` of `lsst.afw.table.ReferenceMatch`
221 List of sources matched to references.
222 wcs : `lsst.afw.geom.SkyWcs`
227 sipObject : `lsst.meas.astrom.sip.CreateWcsWithSip`
230 for i
in range(self.config.numIter):
232 wcs = sipObject.getNewWcs()
236 """Flag deviant matches
238 We return a boolean numpy array indicating whether the corresponding
239 match should be rejected. The previous list of rejections is used
240 so we can calculate uncontaminated statistics.
244 matches : `list` of `lsst.afw.table.ReferenceMatch`
245 List of sources matched to references.
246 wcs : `lsst.afw.geom.SkyWcs`
248 rejected : array-like of `bool`
249 Array of matches rejected from the fit. Unused.
253 rejectedMatches : `ndarray` of type `bool`
254 Matched objects found to be outside of tolerance.
256 fit = [wcs.skyToPixel(m.first.getCoord())
for m
in matches]
257 dx = np.array([ff.getX() - mm.second.getCentroid().getX()
for ff, mm
in zip(fit, matches)])
258 dy = np.array([ff.getY() - mm.second.getCentroid().getY()
for ff, mm
in zip(fit, matches)])
259 good = np.logical_not(rejected)
260 return (dx > self.config.rejSigma*dx[good].
std()) | (dy > self.config.rejSigma*dy[good].
std())
265 We create four plots, for all combinations of (dx, dy) against
266 (x, y). Good points are black, while rejected points are red.
270 matches : `list` of `lsst.afw.table.ReferenceMatch`
271 List of sources matched to references.
272 wcs : `lsst.afw.geom.SkyWcs`
274 rejected : array-like of `bool`
275 Array of matches rejected from the fit.
278 import matplotlib.pyplot
as plt
279 except ImportError
as e:
280 self.log.
warning(
"Unable to import matplotlib: %s", e)
283 fit = [wcs.skyToPixel(m.first.getCoord())
for m
in matches]
284 x1 = np.array([ff.getX()
for ff
in fit])
285 y1 = np.array([ff.getY()
for ff
in fit])
286 x2 = np.array([m.second.getCentroid().getX()
for m
in matches])
287 y2 = np.array([m.second.getCentroid().getY()
for m
in matches])
292 good = np.logical_not(rejected)
294 figure = plt.figure()
295 axes = figure.add_subplot(2, 2, 1)
296 axes.plot(x2[good], dx[good],
'ko')
297 axes.plot(x2[rejected], dx[rejected],
'ro')
299 axes.set_ylabel(
"dx")
301 axes = figure.add_subplot(2, 2, 2)
302 axes.plot(x2[good], dy[good],
'ko')
303 axes.plot(x2[rejected], dy[rejected],
'ro')
305 axes.set_ylabel(
"dy")
307 axes = figure.add_subplot(2, 2, 3)
308 axes.plot(y2[good], dx[good],
'ko')
309 axes.plot(y2[rejected], dx[rejected],
'ro')
311 axes.set_ylabel(
"dx")
313 axes = figure.add_subplot(2, 2, 4)
314 axes.plot(y2[good], dy[good],
'ko')
315 axes.plot(y2[rejected], dy[rejected],
'ro')
317 axes.set_ylabel(
"dy")
An integer coordinate rectangle.
Point in an unspecified spherical coordinate system.
def _fitWcs(self, matches, wcs)
def initialWcs(self, matches, wcs)
def fitWcs(self, matches, initWcs, bbox=None, refCat=None, sourceCat=None, exposure=None)
def plotFit(self, matches, wcs, rejected)
def rejectMatches(self, matches, wcs, rejected)
Vector3d is a vector in ℝ³ with components stored in double precision.
std::shared_ptr< SkyWcs > makeSkyWcs(daf::base::PropertySet &metadata, bool strip=false)
Construct a SkyWcs from FITS keywords.
void updateRefCentroids(geom::SkyWcs const &wcs, ReferenceCollection &refList)
Update centroids in a collection of reference objects.
void updateSourceCoords(geom::SkyWcs const &wcs, SourceCollection &sourceList)
Update sky coordinates in a collection of source objects.
def setMatchDistance(matches)
CreateWcsWithSip< MatchT > makeCreateWcsWithSip(std::vector< MatchT > const &matches, afw::geom::SkyWcs const &linearWcs, int const order, geom::Box2I const &bbox=geom::Box2I(), int const ngrid=0)
Factory function for CreateWcsWithSip.