22 __all__ = [
"ANetAstrometryConfig",
"ANetAstrometryTask",
"showAstrometry"]
29 from lsst.afw.table import Point2DKey, CovarianceMatrix2fKey, updateSourceCoords
31 import lsst.pex.config
as pexConfig
35 from .anetBasicAstrometry
import ANetBasicAstrometryTask
39 solver = pexConfig.ConfigurableField(
40 target=ANetBasicAstrometryTask,
41 doc=
"Basic astrometry solver",
43 forceKnownWcs = pexConfig.Field(dtype=bool, doc=(
44 "Assume that the input image's WCS is correct, without comparing it to any external reality." 45 " (In contrast to using Astrometry.net). NOTE, if you set this, you probably also want to" 46 " un-set 'solver.calculateSip'; otherwise we'll still try to find a TAN-SIP WCS starting " 47 " from the existing WCS"), default=
False)
48 rejectThresh = pexConfig.RangeField(dtype=float, default=3.0, doc=
"Rejection threshold for Wcs fitting",
49 min=0.0, inclusiveMin=
False)
50 rejectIter = pexConfig.RangeField(dtype=int, default=3, doc=
"Rejection iterations for Wcs fitting",
55 """An alias, for a uniform interface with the standard AstrometryTask""" 67 r"""!Use astrometry.net to match input sources with a reference catalog and solve for the Wcs 69 @anchor ANetAstrometryTask_ 71 The actual matching and solving is done by the 'solver'; this Task 72 serves as a wrapper for taking into account the known optical distortion. 74 \section pipe_tasks_astrometry_Contents Contents 76 - \ref pipe_tasks_astrometry_Purpose 77 - \ref pipe_tasks_astrometry_Initialize 78 - \ref pipe_tasks_astrometry_IO 79 - \ref pipe_tasks_astrometry_Config 80 - \ref pipe_tasks_astrometry_Debug 81 - \ref pipe_tasks_astrometry_Example 83 \section pipe_tasks_astrometry_Purpose Description 85 \copybrief ANetAstrometryTask 87 \section pipe_tasks_astrometry_Initialize Task initialisation 91 \section pipe_tasks_astrometry_IO Invoking the Task 95 \section pipe_tasks_astrometry_Config Configuration parameters 97 See \ref ANetAstrometryConfig 99 \section pipe_tasks_astrometry_Debug Debug variables 101 The \link lsst.pipe.base.cmdLineTask.CmdLineTask command line task\endlink interface supports a 102 flag \c -d to import \b debug.py from your \c PYTHONPATH; 103 see \ref baseDebug for more about \b debug.py files. 105 The available variables in ANetAstrometryTask are: 108 <DD> If True call showAstrometry while iterating ANetAstrometryConfig.rejectIter times, 109 and also after converging; and call displayAstrometry after applying the distortion correction. 111 <DD> display frame to use in showAstrometry and displayAstrometry 113 <DD> Pause after showAstrometry and displayAstrometry? 116 \section pipe_tasks_astrometry_Example A complete example of using ANetAstrometryTask 118 See \ref pipe_tasks_photocal_Example. 120 To investigate the \ref pipe_tasks_astrometry_Debug, put something like 124 di = lsstDebug.getInfo(name) # N.b. lsstDebug.Info(name) would call us recursively 125 if name in ("lsst.pipe.tasks.anetAstrometry", "lsst.pipe.tasks.anetBasicAstrometry"): 132 lsstDebug.Info = DebugInfo 134 into your debug.py file and run photoCalTask.py with the \c --debug flag. 136 ConfigClass = ANetAstrometryConfig
137 _DefaultName =
"astrometricSolver" 139 def __init__(self, schema, refObjLoader=None, **kwds):
140 r"""!Create the astrometric calibration task. Most arguments are simply passed onto pipe.base.Task. 142 \param schema An lsst::afw::table::Schema used to create the output lsst.afw.table.SourceCatalog 143 \param refObjLoader The AstrometryTask constructor requires a refObjLoader. In order to make this 144 task retargettable for AstrometryTask it needs to take the same arguments. This argument will be 145 ignored since it uses its own internal loader. 146 \param **kwds keyword arguments to be passed to the lsst.pipe.base.task.Task constructor 148 A centroid field "centroid.distorted" (used internally during the Task's operation) 149 will be added to the schema. 151 pipeBase.Task.__init__(self, **kwds)
154 doc=
"centroid distorted for astrometry solver")
156 doc=
"centroid distorted for astrometry solver")
158 doc=
"centroid distorted err for astrometry solver")
160 doc=
"centroid distorted err for astrometry solver")
162 doc=
"centroid distorted flag astrometry solver")
169 def run(self, exposure, sourceCat):
170 """!Load reference objects, match sources and optionally fit a WCS 172 This is a thin layer around solve or loadAndMatch, depending on config.forceKnownWcs 174 @param[in,out] exposure exposure whose WCS is to be fit 175 The following are read only: 177 - photoCalib (may be absent) 178 - filter (may be unset) 179 - detector (if wcs is pure tangent; may be absent) 180 The following are updated: 181 - wcs (the initial value is used as an initial guess, and is required) 182 @param[in] sourceCat catalog of sourceCat detected on the exposure (an lsst.afw.table.SourceCatalog) 183 @return an lsst.pipe.base.Struct with these fields: 184 - refCat reference object catalog of objects that overlap the exposure (with some margin) 185 (an lsst::afw::table::SimpleCatalog) 186 - matches astrometric matches, a list of lsst.afw.table.ReferenceMatch 187 - matchMeta metadata about the field (an lsst.daf.base.PropertyList) 189 if self.config.forceKnownWcs:
190 return self.
loadAndMatch(exposure=exposure, sourceCat=sourceCat)
192 return self.
solve(exposure=exposure, sourceCat=sourceCat)
195 def solve(self, exposure, sourceCat):
196 r"""!Match with reference sources and calculate an astrometric solution 198 \param[in,out] exposure Exposure to calibrate; wcs is updated 199 \param[in] sourceCat catalog of measured sources (an lsst.afw.table.SourceCatalog) 200 \return a pipeBase.Struct with fields: 201 - refCat reference object catalog of objects that overlap the exposure (with some margin) 202 (an lsst::afw::table::SimpleCatalog) 203 - matches astrometric matches, a list of lsst.afw.table.ReferenceMatch 204 - matchMeta metadata about the field (an lsst.daf.base.PropertyList) 206 The reference catalog actually used is up to the implementation 207 of the solver; it will be manifested in the returned matches as 208 a list of lsst.afw.table.ReferenceMatch objects (\em i.e. of lsst.afw.table.Match with 209 \c first being of type lsst.afw.table.SimpleRecord and \c second type lsst.afw.table.SourceRecord --- 210 the reference object and matched object respectively). 213 The input sources have the centroid slot moved to a new column "centroid.distorted" 214 which has the positions corrected for any known optical distortion; 215 the 'solver' (which is instantiated in the 'astrometry' member) 216 should therefore simply use the centroids provided by calling 217 afw.table.Source.getCentroid() on the individual source records. This column \em must 218 be present in the sources table. 220 \note ignores config.forceKnownWcs 222 results = self.
_astrometry(sourceCat=sourceCat, exposure=exposure)
225 self.
refitWcs(sourceCat=sourceCat, exposure=exposure, matches=results.matches)
231 r"""!Calculate distorted source positions 233 CCD images are often affected by optical distortion that makes 234 the astrometric solution higher order than linear. Unfortunately, 235 most (all?) matching algorithms require that the distortion be 236 small or zero, and so it must be removed. We do this by calculating 237 (un-)distorted positions, based on a known optical distortion model 240 The distortion correction moves sources, so we return the distorted bounding box. 242 \param[in] exposure Exposure to process 243 \param[in,out] sourceCat SourceCatalog; getX() and getY() will be used as inputs, 244 with distorted points in "centroid.distorted" field. 245 \return bounding box of distorted exposure 247 detector = exposure.getDetector()
248 pixToTanXYTransform =
None 250 self.log.
warn(
"No detector associated with exposure; assuming null distortion")
252 pixToTanXYTransform = detector.getTransform(PIXELS, TAN_PIXELS)
254 if pixToTanXYTransform
is None:
255 self.log.
info(
"Null distortion correction")
260 return exposure.getBBox()
263 self.log.
info(
"Applying distortion correction")
265 centroid = pixToTanXYTransform.forwardTransform(s.getCentroid())
272 for corner
in detector.getCorners(TAN_PIXELS):
273 bboxD.include(corner)
279 exposure=exposure, frame=frame, pause=pause)
285 """!Load reference objects overlapping an exposure and match to sources detected on that exposure 287 @param[in] exposure exposure whose WCS is to be fit 288 @param[in] sourceCat catalog of sourceCat detected on the exposure (an lsst.afw.table.SourceCatalog) 289 @param[in] bbox bounding box go use for finding reference objects; if None, use exposure's bbox 291 @return an lsst.pipe.base.Struct with these fields: 292 - refCat reference object catalog of objects that overlap the exposure (with some margin) 293 (an lsst::afw::table::SimpleCatalog) 294 - matches astrometric matches, a list of lsst.afw.table.ReferenceMatch 295 - matchMeta metadata about the field (an lsst.daf.base.PropertyList) 297 @note ignores config.forceKnownWcs 299 bbox = exposure.getBBox()
301 self.makeSubtask(
"solver")
303 astrom = self.
solver.useKnownWcs(
310 if astrom
is None or astrom.getWcs()
is None:
311 raise RuntimeError(
"Unable to solve astrometry")
313 matches = astrom.getMatches()
314 matchMeta = astrom.getMatchMetadata()
315 if matches
is None or len(matches) == 0:
316 raise RuntimeError(
"No astrometric matches")
317 self.log.
info(
"%d astrometric matches", len(matches))
322 frame=frame, pause=
False)
324 return pipeBase.Struct(
325 refCat=astrom.refCat,
331 def _astrometry(self, sourceCat, exposure, bbox=None):
332 r"""!Solve astrometry to produce WCS 334 \param[in] sourceCat Sources on exposure, an lsst.afw.table.SourceCatalog 335 \param[in,out] exposure Exposure to process, an lsst.afw.image.ExposureF or D; wcs is updated 336 \param[in] bbox Bounding box, or None to use exposure 337 \return a pipe.base.Struct with fields: 338 - refCat reference object catalog of objects that overlap the exposure (with some margin) 339 (an lsst::afw::table::SimpleCatalog) 340 - matches astrometric matches, a list of lsst.afw.table.ReferenceMatch 341 - matchMeta metadata about the field (an lsst.daf.base.PropertyList) 343 self.log.
info(
"Solving astrometry")
345 bbox = exposure.getBBox()
348 self.makeSubtask(
"solver")
350 astrom = self.
solver.determineWcs(sourceCat=sourceCat, exposure=exposure, bbox=bbox)
352 if astrom
is None or astrom.getWcs()
is None:
353 raise RuntimeError(
"Unable to solve astrometry")
355 matches = astrom.getMatches()
356 matchMeta = astrom.getMatchMetadata()
357 if matches
is None or len(matches) == 0:
358 raise RuntimeError(
"No astrometric matches")
359 self.log.
info(
"%d astrometric matches", len(matches))
362 exposure.setWcs(astrom.getWcs())
367 frame=frame, pause=
False)
369 return pipeBase.Struct(
370 refCat=astrom.refCat,
377 """!A final Wcs solution after matching and removing distortion 379 Specifically, fitting the non-linear part, since the linear 380 part has been provided by the matching engine. 382 @param sourceCat Sources on exposure, an lsst.afw.table.SourceCatalog 383 @param exposure Exposure of interest, an lsst.afw.image.ExposureF or D 384 @param matches Astrometric matches, as a list of lsst.afw.table.ReferenceMatch 386 @return the resolved-Wcs object, or None if config.solver.calculateSip is False. 389 if self.config.solver.calculateSip:
390 self.log.
info(
"Refitting WCS")
391 origMatches = matches
392 wcs = exposure.getWcs()
399 def fitWcs(initialWcs, title=None):
400 """!Do the WCS fitting and display of the results""" 402 resultWcs = sip.getNewWcs()
404 showAstrometry(exposure, resultWcs, origMatches, matches, frame=frame,
405 title=title, pause=pause)
406 return resultWcs, sip.getScatterOnSky()
410 for i
in range(self.config.rejectIter):
411 wcs, scatter = fitWcs(wcs, title=
"Iteration %d" % i)
413 ref = np.array([wcs.skyToPixel(m.first.getCoord())
for m
in matches])
414 src = np.array([m.second.getCentroid()
for m
in matches])
418 for d, m
in zip(diff, matches):
419 if np.all(np.abs(d) < self.config.rejectThresh*rms):
423 if len(matches) == len(trimmed):
428 wcs, scatter = fitWcs(wcs, title=
"Final astrometry")
431 self.log.
warn(
"Unable to fit SIP: %s", e)
433 self.log.
info(
"Astrometric scatter: %f arcsec (%d matches, %d rejected)",
434 scatter.asArcseconds(), len(matches), numRejected)
440 self.log.
warn(
"Not calculating a SIP solution; matches may be suspect")
445 frame=frame, pause=
False)
450 def showAstrometry(exposure, wcs, allMatches, useMatches, frame=0, title=None, pause=False):
451 r"""!Show results of astrometry fitting 453 \param exposure Image to display 454 \param wcs Astrometric solution 455 \param allMatches List of all astrometric matches (including rejects) 456 \param useMatches List of used astrometric matches 457 \param frame Frame number for display 458 \param title Title for display 459 \param pause Pause to allow viewing of the display and optional debugging? 461 - Matches are shown in yellow if used in the Wcs solution, otherwise red 462 - +: Detected objects 463 - x: Catalogue objects 466 disp = afwDisplay.Display(frame=frame)
467 disp.mtv(exposure, title=title)
469 useIndices =
set(m.second.getId()
for m
in useMatches)
472 with disp.Buffering():
473 for i, m
in enumerate(allMatches):
474 x, y = m.second.getX(), m.second.getY()
475 pix = wcs.skyToPixel(m.first.getCoord())
477 isUsed = m.second.getId()
in useIndices
479 radii.append(np.hypot(pix[0] - x, pix[1] - y))
481 color = afwDisplay.YELLOW
if isUsed
else afwDisplay.RED
483 disp.dot(
"+", x, y, size=10, ctype=color)
484 disp.dot(
"x", pix[0], pix[1], size=10, ctype=color)
486 radii = np.array(radii)
487 print(
"<dr> = %.4g +- %.4g pixels [%d/%d matches]" % (radii.mean(), radii.std(),
488 len(useMatches), len(allMatches)))
494 reply = input(
"Debugging? [p]db [q]uit; any other key to continue... ").
strip()
A floating-point coordinate rectangle geometry.
Reports attempts to exceed implementation-defined length limits for some classes. ...
def loadAndMatch(self, exposure, sourceCat, bbox=None)
Load reference objects overlapping an exposure and match to sources detected on that exposure...
def solve(self, exposure, sourceCat)
Match with reference sources and calculate an astrometric solution.
PointKey< double > Point2DKey
def distort(self, sourceCat, exposure)
Calculate distorted source positions.
def showAstrometry(exposure, wcs, allMatches, useMatches, frame=0, title=None, pause=False)
Show results of astrometry fitting.
daf::base::PropertySet * set
void updateSourceCoords(geom::SkyWcs const &wcs, SourceCollection &sourceList)
Update sky coordinates in a collection of source objects.
def refitWcs(self, sourceCat, exposure, matches)
A final Wcs solution after matching and removing distortion.
def _astrometry(self, sourceCat, exposure, bbox=None)
Solve astrometry to produce WCS.
def run(self, exposure, sourceCat)
Load reference objects, match sources and optionally fit a WCS.
def displayAstrometry(refCat=None, sourceCat=None, distortedCentroidKey=None, bbox=None, exposure=None, matches=None, frame=1, title="", pause=True)
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.
An integer coordinate rectangle.
def __init__(self, schema, refObjLoader=None, kwds)
Create the astrometric calibration task.
Use astrometry.net to match input sources with a reference catalog and solve for the Wcs...