diff --git a/python/lsst/ip/diffim/getTemplate.py b/python/lsst/ip/diffim/getTemplate.py index 70382b03..6d03a25b 100644 --- a/python/lsst/ip/diffim/getTemplate.py +++ b/python/lsst/ip/diffim/getTemplate.py @@ -243,6 +243,7 @@ def getExposures(self, coaddExposureHandles, bbox, skymap, wcs): # Exposure's validPolygon would be more accurate detectorPolygon = geom.Box2D(bbox) + detectorCorners = wcs.pixelToSky(detectorPolygon.getCorners()) overlappingArea = 0 coaddExposures = collections.defaultdict(list) dataIds = collections.defaultdict(list) @@ -251,11 +252,15 @@ def getExposures(self, coaddExposureHandles, bbox, skymap, wcs): dataId = coaddRef.dataId patchWcs = skymap[dataId["tract"]].getWcs() patchBBox = skymap[dataId["tract"]][dataId["patch"]].getOuterBBox() - patchCorners = patchWcs.pixelToSky(geom.Box2D(patchBBox).getCorners()) - patchPolygon = afwGeom.Polygon(wcs.skyToPixel(patchCorners)) - if patchPolygon.intersection(detectorPolygon): + patchPolygon = afwGeom.Polygon(geom.Box2D(patchBBox)) + # Calculate detector/patch overlap in patch coordinates rather than + # detector coordinates because the skymap's inverse mapping + # (patchWcs.skyToPixel()) is more stable than the detector's for + # arbitrary sky coordinates. + detectorInPatchCoordinates = afwGeom.Polygon(patchWcs.skyToPixel(detectorCorners)) + if patchPolygon.intersection(detectorInPatchCoordinates): overlappingArea += patchPolygon.intersectionSingle( - detectorPolygon + detectorInPatchCoordinates ).calculateArea() self.log.info( "Using template input tract=%s, patch=%s",