Merge pull request #102 from lsst/tickets/DM-15635

DM-15635: Approximate filter throughput for DcrCoadds

python/lsst/ip/diffim/dcrModel.py

@@ -289,7 +289,8 @@ def assign(self, dcrSubModel, bbox=None):
             model.assign(subModel[bbox], bbox)
 
     def buildMatchedTemplate(self, exposure=None, warpCtrl=None,
-                             visitInfo=None, bbox=None, wcs=None, mask=None):
+                             visitInfo=None, bbox=None, wcs=None, mask=None,
+                             splitSubfilters=False):
         """Create a DCR-matched template image for an exposure.
 
         Parameters
@@ -310,6 +311,9 @@ def buildMatchedTemplate(self, exposure=None, warpCtrl=None,
             Ignored if ``exposure`` is set.
         mask : `lsst.afw.image.Mask`, optional
             reference mask to use for the template image.
+        splitSubfilters : `bool`, optional
+            Calculate DCR for two evenly-spaced wavelengths in each subfilter,
+            instead of at the midpoint. Default: False
 
         Returns
         -------
@@ -335,10 +339,10 @@ def buildMatchedTemplate(self, exposure=None, warpCtrl=None,
             warpCtrl = afwMath.WarpingControl("lanczos3", "bilinear",
                                               cacheSize=0, interpLength=max(bbox.getDimensions()))
 
-        dcrShift = calculateDcr(visitInfo, wcs, self.filter, len(self))
+        dcrShift = calculateDcr(visitInfo, wcs, self.filter, len(self), splitSubfilters=splitSubfilters)
         templateImage = afwImage.MaskedImageF(bbox)
         for subfilter, dcr in enumerate(dcrShift):
-            templateImage += applyDcr(self[subfilter][bbox], dcr, warpCtrl)
+            templateImage += applyDcr(self[subfilter][bbox], dcr, warpCtrl, splitSubfilters=splitSubfilters)
         if mask is not None:
             templateImage.setMask(mask[bbox])
         return templateImage
@@ -533,7 +537,7 @@ def applyImageThresholds(self, maskedImage, highThreshold=None, lowThreshold=Non
             image[lowPixels] = lowThreshold[lowPixels]
 
 
-def applyDcr(maskedImage, dcr, warpCtrl, bbox=None, useInverse=False):
+def applyDcr(maskedImage, dcr, warpCtrl, bbox=None, useInverse=False, splitSubfilters=False):
     """Shift a masked image.
 
     Parameters
@@ -549,23 +553,38 @@ def applyDcr(maskedImage, dcr, warpCtrl, bbox=None, useInverse=False):
         Shifts the entire image if None (Default).
     useInverse : `bool`, optional
         Use the reverse of ``dcr`` for the shift. Default: False
+    splitSubfilters : `bool`, optional
+        Calculate DCR for two evenly-spaced wavelengths in each subfilter,
+        instead of at the midpoint. Default: False
 
     Returns
     -------
-    `lsst.afw.image.maskedImageF`
+    shiftedImage : `lsst.afw.image.maskedImageF`
         A masked image, with the pixels within the bounding box shifted.
     """
     padValue = afwImage.pixel.SinglePixelF(0., maskedImage.mask.getPlaneBitMask("NO_DATA"), 0)
     if bbox is None:
         bbox = maskedImage.getBBox()
-    shiftedImage = afwImage.MaskedImageF(bbox)
-    transform = makeTransform(AffineTransform((-1.0 if useInverse else 1.0)*dcr))
-    afwMath.warpImage(shiftedImage, maskedImage[bbox],
-                      transform, warpCtrl, padValue=padValue)
+    if splitSubfilters:
+        shiftedImage = afwImage.MaskedImageF(bbox)
+        transform0 = makeTransform(AffineTransform((-1.0 if useInverse else 1.0)*dcr[0]))
+        afwMath.warpImage(shiftedImage, maskedImage[bbox],
+                          transform0, warpCtrl, padValue=padValue)
+        shiftedImage1 = afwImage.MaskedImageF(bbox)
+        transform1 = makeTransform(AffineTransform((-1.0 if useInverse else 1.0)*dcr[1]))
+        afwMath.warpImage(shiftedImage1, maskedImage[bbox],
+                          transform1, warpCtrl, padValue=padValue)
+        shiftedImage += shiftedImage1
+        shiftedImage /= 2.
+    else:
+        shiftedImage = afwImage.MaskedImageF(bbox)
+        transform = makeTransform(AffineTransform((-1.0 if useInverse else 1.0)*dcr))
+        afwMath.warpImage(shiftedImage, maskedImage[bbox],
+                          transform, warpCtrl, padValue=padValue)
     return shiftedImage
 
 
-def calculateDcr(visitInfo, wcs, filterInfo, dcrNumSubfilters):
+def calculateDcr(visitInfo, wcs, filterInfo, dcrNumSubfilters, splitSubfilters=False):
     """Calculate the shift in pixels of an exposure due to DCR.
 
     Parameters
@@ -578,14 +597,18 @@ def calculateDcr(visitInfo, wcs, filterInfo, dcrNumSubfilters):
         The filter definition, set in the current instruments' obs package.
     dcrNumSubfilters : `int`
         Number of sub-filters used to model chromatic effects within a band.
+    splitSubfilters : `bool`, optional
+        Calculate DCR for two evenly-spaced wavelengths in each subfilter,
+        instead of at the midpoint. Default: False
 
     Returns
     -------
-    `lsst.afw.geom.Extent2I`
+    dcrShift : `lsst.afw.geom.Extent2I`
         The 2D shift due to DCR, in pixels.
     """
     rotation = calculateImageParallacticAngle(visitInfo, wcs)
     dcrShift = []
+    weight = [0.75, 0.25]
     lambdaEff = filterInfo.getFilterProperty().getLambdaEff()
     for wl0, wl1 in wavelengthGenerator(filterInfo, dcrNumSubfilters):
         # Note that diffRefractAmp can be negative, since it's relative to the midpoint of the full band
@@ -597,11 +620,20 @@ def calculateDcr(visitInfo, wcs, filterInfo, dcrNumSubfilters):
                                                  elevation=visitInfo.getBoresightAzAlt().getLatitude(),
                                                  observatory=visitInfo.getObservatory(),
                                                  weather=visitInfo.getWeather())
-        diffRefractAmp = (diffRefractAmp0 + diffRefractAmp1)/2.
-        diffRefractPix = diffRefractAmp.asArcseconds()/wcs.getPixelScale().asArcseconds()
-        shiftX = diffRefractPix*np.sin(rotation.asRadians())
-        shiftY = diffRefractPix*np.cos(rotation.asRadians())
-        dcrShift.append(afwGeom.Extent2D(shiftX, shiftY))
+        if splitSubfilters:
+            diffRefractPix0 = diffRefractAmp0.asArcseconds()/wcs.getPixelScale().asArcseconds()
+            diffRefractPix1 = diffRefractAmp1.asArcseconds()/wcs.getPixelScale().asArcseconds()
+            diffRefractArr = [diffRefractPix0*weight[0] + diffRefractPix1*weight[1],
+                              diffRefractPix0*weight[1] + diffRefractPix1*weight[0]]
+            shiftX = [diffRefractPix*np.sin(rotation.asRadians()) for diffRefractPix in diffRefractArr]
+            shiftY = [diffRefractPix*np.cos(rotation.asRadians()) for diffRefractPix in diffRefractArr]
+            dcrShift.append((afwGeom.Extent2D(shiftX[0], shiftY[0]), afwGeom.Extent2D(shiftX[1], shiftY[1])))
+        else:
+            diffRefractAmp = (diffRefractAmp0 + diffRefractAmp1)/2.
+            diffRefractPix = diffRefractAmp.asArcseconds()/wcs.getPixelScale().asArcseconds()
+            shiftX = diffRefractPix*np.sin(rotation.asRadians())
+            shiftY = diffRefractPix*np.cos(rotation.asRadians())
+            dcrShift.append(afwGeom.Extent2D(shiftX, shiftY))
     return dcrShift