Merge pull request #221 from lsst/tickets/DM-36768

DM-36768: Add psf image to alert cutouts

python/lsst/ap/association/packageAlerts.py

@@ -265,23 +265,28 @@ def run(self,
  sphPoint = geom.SpherePoint(diaSource["ra"],
  diaSource["dec"],
  geom.degrees)
+ pixelPoint = geom.Point2D(diaSource["x"],
+ diaSource["y"])
 
  cutoutExtent = self.createDiaSourceExtent(diaSource["bboxSize"])
  diffImCutout = self.createCcdDataCutout(
  diffIm,
  sphPoint,
+ pixelPoint,
  cutoutExtent,
  diffImPhotoCalib,
  diaSource["diaSourceId"])
  calexpCutout = self.createCcdDataCutout(
  calexp,
  sphPoint,
+ pixelPoint,
  cutoutExtent,
  calexpPhotoCalib,
  diaSource["diaSourceId"])
  templateCutout = self.createCcdDataCutout(
  template,
  sphPoint,
+ pixelPoint,
  cutoutExtent,
  templatePhotoCalib,
  diaSource["diaSourceId"])
@@ -370,7 +375,7 @@ def produceAlerts(self, alerts, visit, detector):
 
  self.producer.flush()
 
- def createCcdDataCutout(self, image, skyCenter, extent, photoCalib, srcId):
+ def createCcdDataCutout(self, image, skyCenter, pixelCenter, extent, photoCalib, srcId):
  """Grab an image as a cutout and return a calibrated CCDData image.
 
  Parameters
@@ -379,6 +384,8 @@ def createCcdDataCutout(self, image, skyCenter, extent, photoCalib, srcId):
  Image to pull cutout from.
  skyCenter : `lsst.geom.SpherePoint`
  Center point of DiaSource on the sky.
+ pixelCenter : `lsst.geom.Point2D`
+ Pixel center of DiaSource on the sky.
  extent : `lsst.geom.Extent2I`
  Bounding box to cutout from the image.
  photoCalib : `lsst.afw.image.PhotoCalib`
@@ -393,18 +400,18 @@ def createCcdDataCutout(self, image, skyCenter, extent, photoCalib, srcId):
  CCDData object storing the calibrate information from the input
  difference or template image.
  """
+
  # Catch errors in retrieving the cutout.
  try:
- cutout = image.getCutout(skyCenter, extent)
+ cutout = image.getCutout(pixelCenter, extent)
  except InvalidParameterError:
- point = image.getWcs().skyToPixel(skyCenter)
  imBBox = image.getBBox()
- if not geom.Box2D(image.getBBox()).contains(point):
+ if not geom.Box2D(image.getBBox()).contains(pixelCenter):
  self.log.warning(
  "DiaSource id=%i centroid lies at pixel (%.2f, %.2f) "
  "which is outside the Exposure with bounding box "
  "((%i, %i), (%i, %i)). Returning None for cutout...",
- srcId, point.x, point.y,
+ srcId, pixelCenter.x, pixelCenter.y,
  imBBox.minX, imBBox.maxX, imBBox.minY, imBBox.maxY)
  else:
  raise InvalidParameterError(
@@ -414,6 +421,9 @@ def createCcdDataCutout(self, image, skyCenter, extent, photoCalib, srcId):
  % srcId)
  return None
 
+ # use image.psf.computeKernelImage to provide PSF centered in the array
+ cutoutPsf = image.psf.computeKernelImage(pixelCenter).array
+
  # Find the value of the bottom corner of our cutout's BBox and
  # subtract 1 so that the CCDData cutout position value will be
  # [1, 1].
@@ -437,6 +447,7 @@ def createCcdDataCutout(self, image, skyCenter, extent, photoCalib, srcId):
  uncertainty=VarianceUncertainty(calibCutout.getVariance().array),
  flags=calibCutout.getMask().array,
  wcs=cutoutWcs,
+ psf=cutoutPsf,
  meta={"cutMinX": cutOutMinX,
  "cutMinY": cutOutMinY},
  unit=u.nJy)

tests/test_packageAlerts.py

@@ -255,6 +255,7 @@ def testCreateCcdDataCutout(self):
  ccdData = packageAlerts.createCcdDataCutout(
  self.exposure,
  self.exposure.getWcs().getSkyOrigin(),
+ self.exposure.getWcs().getPixelOrigin(),
  self.exposure.getBBox().getDimensions(),
  self.exposure.getPhotoCalib(),
  diaSrcId)
@@ -265,10 +266,13 @@ def testCreateCcdDataCutout(self):
  self.cutoutWcs.wcs.cd)
  self.assertFloatsAlmostEqual(ccdData.data,
  calibExposure.getImage().array)
+ self.assertFloatsAlmostEqual(ccdData.psf,
+ self.exposure.psf.computeKernelImage(self.center).array)
 
  ccdData = packageAlerts.createCcdDataCutout(
  self.exposure,
  geom.SpherePoint(0, 0, geom.degrees),
+ geom.Point2D(-100, -100),
  self.exposure.getBBox().getDimensions(),
  self.exposure.getPhotoCalib(),
  diaSrcId)
@@ -321,12 +325,14 @@ def testMakeAlertDict(self):
  sphPoint = geom.SpherePoint(diaSource["ra"],
  diaSource["dec"],
  geom.degrees)
+ pixelPoint = geom.Point2D(diaSource["x"], diaSource["y"])
  cutout = self.exposure.getCutout(sphPoint,
  geom.Extent2I(self.cutoutSize,
  self.cutoutSize))
  ccdCutout = packageAlerts.createCcdDataCutout(
  cutout,
  sphPoint,
+ pixelPoint,
  geom.Extent2I(self.cutoutSize, self.cutoutSize),
  cutout.getPhotoCalib(),
  1234)
@@ -509,12 +515,14 @@ def testRun_without_produce(self, mock_server_check):
  sphPoint = geom.SpherePoint(alert["diaSource"]["ra"],
  alert["diaSource"]["dec"],
  geom.degrees)
+ pixelPoint = geom.Point2D(alert["diaSource"]["x"], alert["diaSource"]["y"])
  cutout = self.exposure.getCutout(sphPoint,
  geom.Extent2I(self.cutoutSize,
  self.cutoutSize))
  ccdCutout = packageAlerts.createCcdDataCutout(
  cutout,
  sphPoint,
+ pixelPoint,
  geom.Extent2I(self.cutoutSize, self.cutoutSize),
  cutout.getPhotoCalib(),
  1234)