Skip to content

Commit

Permalink
Add maxMeanOffsetArcsec config to impose WCS fail
Browse files Browse the repository at this point in the history
It is desirable to have a threshold above which even a converged WCS
fit will be deemed a failure, resulting in the WCS attached to the
exposure being set to None and all coord_ra and coord_dec entries set
to np.nan in the catalogs.  Downstream processing can then appropriately
decide whether this exposure is suitable for further consideration (e.g.
it will be skipped in the coaddition process).

The specific value for this threshold will depend on the data and
processing workflow.  For example, if the single frame measurement is
the final WCS to be associated with an exposure before moving on to the
coaddition stage, one might want to be fairly strict with this threshold
so as not to allow exposures with "converged" but "bad" WCS fits making
there way into a coadd.  A rough upper bound for inputs to coadds of
about half a pixel has been suggested in this case.  On the other hand,
if the workflow includes an external (and global) astrometric fitting step
(e.g. jointcal) , then it may be possible to loosen this constraint in the
hopes of allowing the external calibration to recover a good WCS fit for a
"bad" initial fit from single frame processing.

Since the main Rubin/LSST processing workflow will indeed include an
external/global calibration step, we set the default here to match the
maximum possible "badness" of the SFM WCS that it could conceivably
recover from, which is currently on the order of 1 arcsec.
  • Loading branch information
laurenam committed Aug 31, 2021
1 parent 2d91323 commit b8397f1
Showing 1 changed file with 61 additions and 44 deletions.
105 changes: 61 additions & 44 deletions python/lsst/meas/astrom/astrometry.py
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,15 @@ class AstrometryConfig(RefMatchConfig):
default=0.001,
min=0,
)
maxMeanOffsetArcsec = pexConfig.RangeField(
doc="Maximum mean offset (in arcsec) of a WCS fit beyond which the fit is considered "
"a failure. Default is set to maximum tolerated by the external global calibration "
"(i.e. jointcal) step for conceivable recovery. Appropriate value will be specific "
"dataset and workflow dependent.",
dtype=float,
default=1.0,
min=0,
)
doMagnitudeOutlierRejection = pexConfig.Field(
dtype=bool,
doc=("If True then a rough zeropoint will be computed from matched sources "
Expand Down Expand Up @@ -238,47 +247,60 @@ def solve(self, exposure, sourceCat):
match_tolerance = None
fitFailed = False
for i in range(self.config.maxIter):
iterNum = i + 1
try:
tryRes = self._matchAndFitWcs(
refCat=refSelection.sourceCat,
sourceCat=sourceCat,
goodSourceCat=sourceSelection.sourceCat,
refFluxField=loadRes.fluxField,
bbox=expMd.bbox,
wcs=wcs,
exposure=exposure,
match_tolerance=match_tolerance,
)
except Exception as e:
# If we have had a succeessful iteration then use that;
# otherwise fail.
if i > 0:
self.log.info("Fit WCS iter %d failed; using previous iteration: %s" % (iterNum, e))
iterNum -= 1
break
else:
fitFailed = True

if not fitFailed:
match_tolerance = tryRes.match_tolerance
tryMatchDist = self._computeMatchStatsOnSky(tryRes.matches)
self.log.debug(
"Match and fit WCS iteration %d: found %d matches with mean and scatter = "
"%0.3f +- %0.3f arcsec; max match distance = %0.3f arcsec",
iterNum, len(tryRes.matches), tryMatchDist.distMean.asArcseconds(),
tryMatchDist.distStdDev.asArcseconds(), tryMatchDist.maxMatchDist.asArcseconds())

maxMatchDist = tryMatchDist.maxMatchDist
res = tryRes
wcs = res.wcs
if maxMatchDist.asArcseconds() < self.config.minMatchDistanceArcSec:
iterNum = i + 1
try:
tryRes = self._matchAndFitWcs(
refCat=refSelection.sourceCat,
sourceCat=sourceCat,
goodSourceCat=sourceSelection.sourceCat,
refFluxField=loadRes.fluxField,
bbox=expMd.bbox,
wcs=wcs,
exposure=exposure,
match_tolerance=match_tolerance,
)
except Exception as e:
# If we have had a succeessful iteration then use that;
# otherwise fail.
if i > 0:
self.log.info("Fit WCS iter %d failed; using previous iteration: %s" % (iterNum, e))
iterNum -= 1
break
else:
self.log.info("Fit WCS iter %d failed: %s" % (iterNum, e))
fitFailed = True

if not fitFailed:
match_tolerance = tryRes.match_tolerance
tryMatchDist = self._computeMatchStatsOnSky(tryRes.matches)
self.log.debug(
"Max match distance = %0.3f arcsec < %0.3f = config.minMatchDistanceArcSec; "
"that's good enough",
maxMatchDist.asArcseconds(), self.config.minMatchDistanceArcSec)
break
match_tolerance.maxMatchDist = maxMatchDist
"Match and fit WCS iteration %d: found %d matches with mean and scatter = "
"%0.3f +- %0.3f arcsec; max match distance = %0.3f arcsec",
iterNum, len(tryRes.matches), tryMatchDist.distMean.asArcseconds(),
tryMatchDist.distStdDev.asArcseconds(), tryMatchDist.maxMatchDist.asArcseconds())

maxMatchDist = tryMatchDist.maxMatchDist
res = tryRes
wcs = res.wcs
if maxMatchDist.asArcseconds() < self.config.minMatchDistanceArcSec:
self.log.debug(
"Max match distance = %0.3f arcsec < %0.3f = config.minMatchDistanceArcSec; "
"that's good enough",
maxMatchDist.asArcseconds(), self.config.minMatchDistanceArcSec)
break
match_tolerance.maxMatchDist = maxMatchDist

if not fitFailed:
self.log.info("Matched and fit WCS in %d iterations; "
"found %d matches with mean and scatter = %0.3f +- %0.3f arcsec" %
(iterNum, len(tryRes.matches), tryMatchDist.distMean.asArcseconds(),
tryMatchDist.distStdDev.asArcseconds()))
if not fitFailed and tryMatchDist.distMean.asArcseconds() > self.config.maxMeanOffsetArcsec:
self.log.info("Assigning as a fit failure: mean offset on sky = %0.3f arcsec > %0.3f "
"config.maxMeanOffsetArcsec" % (tryMatchDist.distMean.asArcseconds(),
self.config.maxMeanOffsetArcsec))
fitFailed = True

if fitFailed:
self.log.warn("WCS fit failed. Setting exposure's WCS to None and coord_ra & coord_dec cols "
Expand All @@ -289,11 +311,6 @@ def solve(self, exposure, sourceCat):
matches = None
scatterOnSky = None
else:
self.log.info(
"Matched and fit WCS in %d iterations; "
"found %d matches with mean and scatter = %0.3f +- %0.3f arcsec" %
(iterNum, len(tryRes.matches), tryMatchDist.distMean.asArcseconds(),
tryMatchDist.distStdDev.asArcseconds()))
for m in res.matches:
if self.usedKey:
m.second.set(self.usedKey, True)
Expand Down

0 comments on commit b8397f1

Please sign in to comment.