Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

DM-24325: Include CCDData cutouts in alert packets #89

Merged
merged 4 commits into from
Jun 22, 2020
Merged

DM-24325: Include CCDData cutouts in alert packets #89

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
18f63c8
Change to CCDData instead of afw in alerts.
morriscb Jun 18, 2020
34fb04a
Add last unittest.
morriscb Jun 19, 2020
865c859
Add/modify doc strings.
morriscb Jun 19, 2020
95ff93e
Respond to review
morriscb Jun 22, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
118 changes: 104 additions & 14 deletions python/lsst/ap/association/packageAlerts.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,10 +21,15 @@

__all__ = ("PackageAlertsConfig", "PackageAlertsTask")

import io
import os

import lsst.afw.fits as afwFits
from astropy import wcs
import astropy.units as u
from astropy.nddata import CCDData, VarianceUncertainty

import lsst.alert.packet as alertPack
import lsst.afw.geom as afwGeom
import lsst.geom as geom
import lsst.pex.config as pexConfig
import lsst.pipe.base as pipeBase
Expand Down Expand Up @@ -66,6 +71,8 @@ class PackageAlertsTask(pipeBase.Task):
ConfigClass = PackageAlertsConfig
_DefaultName = "packageAlerts"

_scale = (1.0 * geom.arcseconds).asDegrees()

def __init__(self, **kwargs):
super().__init__(**kwargs)
self.alertSchema = alertPack.Schema.from_file(self.config.schemaFile)
Expand Down Expand Up @@ -106,6 +113,7 @@ def run(self,
self._patchDiaSources(diaSourceCat)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's not part of your changes, but it doesn't look like template or ccdExposureIdBits are used anywhere.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Eventually will so I'm going to leave it for now.

self._patchDiaSources(diaSrcHistory)
ccdVisitId = diffIm.getInfo().getVisitInfo().getExposureId()
diffImPhotoCalib = diffIm.getPhotoCalib()
for srcIndex, diaSource in diaSourceCat.iterrows():
# Get all diaSources for the associated diaObject.
diaObject = diaObjectCat.loc[srcIndex[0]]
Expand All @@ -117,7 +125,11 @@ def run(self,
diaSource["decl"],
geom.degrees)
cutoutBBox = self.createDiaSourceBBox(diaSource["bboxSize"])
diffImCutout = diffIm.getCutout(sphPoint, cutoutBBox)
diffImCutout = self.createCcdDataCutout(
diffIm.getCutout(sphPoint, cutoutBBox),
sphPoint,
diffImPhotoCalib)

templateCutout = None
# TODO: Create alertIds DM-24858
alertId = diaSource["diaSourceId"]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also not part of your changes, but I don't see how _patchDiaSources below changes the "currently grouped alert data to None" as the documentation says.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removed. Should have been removed in a previous ticket but was not.

Expand All @@ -134,8 +146,7 @@ def run(self,
self.alertSchema.store_alerts(f, alerts)

def _patchDiaSources(self, diaSources):
"""Add the ``programId`` column to the data and change currently
grouped alert data to ``None``.
"""Add the ``programId`` column to the data.

Parameters
----------
Expand Down Expand Up @@ -165,6 +176,84 @@ def createDiaSourceBBox(self, bboxSize):
bbox = geom.Extent2I(bboxSize, bboxSize)
return bbox

def createCcdDataCutout(self, cutout, skyCenter, photoCalib):
"""Convert a cutout into a calibrate CCDData image.

Parameters
----------
cutout : `lsst.afw.image.ExposureF`
Cutout to convert.
skyCenter : `lsst.geom.SpherePoint`
Center point of DiaSource on the sky.
photoCalib : `lsst.afw.image.PhotoCalib`
Calibrate object of the image the cutout is cut from.

Returns
-------
ccdData : `astropy.nddata.CCDData`
CCDData object storing the calibrate information from the input
difference or template image.
"""
# 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].
cutOutMinX = cutout.getBBox().minX - 1
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please add a comment explaining the additional -1 here.

cutOutMinY = cutout.getBBox().minY - 1
center = cutout.getWcs().skyToPixel(skyCenter)
calibCutout = photoCalib.calibrateImage(cutout.getMaskedImage())

cutoutWcs = wcs.WCS(naxis=2)
cutoutWcs.array_shape = (cutout.getBBox().getWidth(),
cutout.getBBox().getWidth())
cutoutWcs.wcs.crpix = [center.x - cutOutMinX, center.y - cutOutMinY]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is the crpix at the corner rather than the center of the cutout?

Copy link
Contributor Author

@morriscb morriscb Jun 22, 2020 •
edited

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So, CCDData doesn't store absolute pixel locations. So when you create a dataset with it the "bottom left" pixel is always [1, 1]. Thus we center the crpix such that it's correctly within this newly shifted cutout and at the location of the DiaSource in the cutout. There is a possible future feature in astropy to add an origin to the CCDData object but it is currently not a priority.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think "not a priority" is right for them but if we want to add that feature they seem like they would be more than happy to accept a patch. We have been in discussion for years to add pixel origin to NDData.

cutoutWcs.wcs.crval = [skyCenter.getRa().asDegrees(),
skyCenter.getDec().asDegrees()]
cutoutWcs.wcs.cd = self.makeLocalTransformMatrix(cutout.getWcs(),
center,
skyCenter)

return CCDData(
data=calibCutout.getImage().array,
uncertainty=VarianceUncertainty(calibCutout.getVariance().array),
flags=calibCutout.getMask().array,
wcs=cutoutWcs,
meta={"cutMinX": cutOutMinX,
"cutMinY": cutOutMinY},
unit=u.nJy)

def makeLocalTransformMatrix(self, wcs, center, skyCenter):
"""Create a local, linear approximation of the wcs transformation
matrix.

The approximation is created as if the center is at RA=0, DEC=0. All
comparing x,y coordinate are relative to the position of center. Matrix
is initially calculated with units arcseconds and then converted to
degrees. This yields higher precision results due to quirks in AST.

Parameters
----------
wcs : `lsst.afw.geom.SkyWcs`
Wcs to approximate
center : `lsst.geom.Point2D`
Point at which to evaluate the LocalWcs.
skyCenter : `lsst.geom.SpherePoint`
Point on sky to approximate the Wcs.

Returns
-------
localMatrix : `numpy.ndarray`
Matrix representation the local wcs approximation with units
degrees.
"""
blankCDMatrix = [[self._scale, 0], [0, self._scale]]
localGnomonicWcs = afwGeom.makeSkyWcs(
center, skyCenter, blankCDMatrix)
measurementToLocalGnomonic = wcs.getTransform().then(
localGnomonicWcs.getTransform().inverted()
)
localMatrix = measurementToLocalGnomonic.getJacobian(center)
return localMatrix / 3600

def makeAlertDict(self,
alertId,
diaSource,
Expand All @@ -182,12 +271,12 @@ def makeAlertDict(self,
DiaObject that ``diaSource`` is matched to.
objDiaSrcHistory : `pandas.DataFrame`
12 month history of ``diaObject`` excluding the latest DiaSource.
diffImCutout : `lsst.afw.image.ExposureF` or `None`
diffImCutout : `astropy.nddata.CCDData` or `None`
Cutout of the difference image around the location of ``diaSource``
with a size set by the ``cutoutSize`` configurable.
templateCutout : `lsst.afw.image.ExposureF` or `None`
with a min size set by the ``cutoutSize`` configurable.
templateCutout : `astropy.nddata.CCDData` or `None`
Cutout of the template image around the location of ``diaSource``
with a size set by the ``cutoutSize`` configurable.
with a min size set by the ``cutoutSize`` configurable.
"""
alert = dict()
alert['alertId'] = alertId
Expand All @@ -205,25 +294,26 @@ def makeAlertDict(self,

alert['ssObject'] = None

alert['cutoutDifference'] = self.makeCutoutBytes(diffImCutout)
alert['cutoutDifference'] = self.streamCcdDataToBytes(diffImCutout)
# TODO: add template cutouts in DM-24327
alert["cutoutTemplate"] = None

return alert

def makeCutoutBytes(self, cutout):
def streamCcdDataToBytes(self, cutout):
"""Serialize a cutout into bytes.

Parameters
----------
cutout : `lsst.afw.image.ExposureF`
cutout : `astropy.nddata.CCDData`
Cutout to serialize.

Returns
-------
coutputBytes : `bytes`
Input cutout serialized into byte data.
"""
mgr = afwFits.MemFileManager()
cutout.writeFits(mgr)
return mgr.getData()
with io.BytesIO() as streamer:
cutout.write(streamer, format="fits")
cutoutBytes = streamer.getvalue()
return cutoutBytes
81 changes: 67 additions & 14 deletions tests/test_packageAlerts.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,19 +19,22 @@
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.

import io
import os
import numpy as np
import pandas as pd
import shutil
import tempfile
import unittest

from astropy import wcs
from astropy.nddata import CCDData

from lsst.ap.association import (PackageAlertsConfig,
PackageAlertsTask,
make_dia_source_schema,
make_dia_object_schema)
from lsst.afw.cameraGeom.testUtils import DetectorWrapper
import lsst.afw.fits as afwFits
import lsst.afw.image as afwImage
import lsst.afw.image.utils as afwImageUtils
import lsst.daf.base as dafBase
Expand Down Expand Up @@ -225,7 +228,7 @@ def _roundTripThroughApdb(objects, sources, dateTime):
return (diaObjects, diaSources)


class TestPackageAlerts(unittest.TestCase):
class TestPackageAlerts(lsst.utils.tests.TestCase):

def setUp(self):
np.random.seed(1234)
Expand Down Expand Up @@ -273,6 +276,14 @@ def setUp(self):
self.diaSourceHistory = diaSourceHistory.drop(labels=[(0, "g", 8),
(1, "g", 9)])

self.cutoutWcs = wcs.WCS(naxis=2)
self.cutoutWcs.wcs.crpix = [self.center[0], self.center[1]]
self.cutoutWcs.wcs.crval = [
self.exposure.getWcs().getSkyOrigin().getRa().asDegrees(),
self.exposure.getWcs().getSkyOrigin().getDec().asDegrees()]
self.cutoutWcs.wcs.cd = self.exposure.getWcs().getCdMatrix()
self.cutoutWcs.wcs.ctype = ["RA---TAN", "DEC--TAN"]

def testCreateBBox(self):
"""Test the bbox creation
"""
Expand All @@ -288,22 +299,59 @@ def testCreateBBox(self):
self.assertTrue(bbox == geom.Extent2I(self.cutoutSize,
self.cutoutSize))

def testMakeCutoutBytes(self):
"""Test round tripping an exposure/cutout to bytes and back.
def testCreateCcdDataCutout(self):
"""Test that the data is being extracted into the CCDData cutout
correctly.
"""
packageAlerts = PackageAlertsTask()

ccdData = packageAlerts.createCcdDataCutout(
self.exposure,
self.exposure.getWcs().getSkyOrigin(),
self.exposure.getPhotoCalib())
calibExposure = self.exposure.getPhotoCalib().calibrateImage(
self.exposure.getMaskedImage())

self.assertFloatsAlmostEqual(ccdData.wcs.wcs.cd,
self.cutoutWcs.wcs.cd)
self.assertFloatsAlmostEqual(ccdData.data,
calibExposure.getImage().array)

def testMakeLocalTransformMatrix(self):
"""Test that the local WCS approximation is correct.
"""
packageAlerts = PackageAlertsTask()

sphPoint = self.exposure.getWcs().pixelToSky(self.center)
cutout = self.exposure.getCutout(sphPoint,
geom.Extent2I(self.cutoutSize,
self.cutoutSize))
cd = packageAlerts.makeLocalTransformMatrix(
cutout.getWcs(), self.center, sphPoint)
self.assertFloatsAlmostEqual(
cd,
cutout.getWcs().getCdMatrix(),
rtol=1e-11,
atol=1e-11)

def testStreamCcdDataToBytes(self):
"""Test round tripping an CCDData cutout to bytes and back.
"""
packageAlerts = PackageAlertsTask()

sphPoint = self.exposure.getWcs().pixelToSky(self.center)
cutout = self.exposure.getCutout(sphPoint,
geom.Extent2I(self.cutoutSize,
self.cutoutSize))
cutoutCcdData = CCDData(
data=cutout.getImage().array,
wcs=self.cutoutWcs,
unit="adu")

cutoutBytes = packageAlerts.makeCutoutBytes(cutout)
tempMemFile = afwFits.MemFileManager(len(cutoutBytes))
tempMemFile.setData(cutoutBytes, len(cutoutBytes))
cutoutFromBytes = afwImage.ExposureF(tempMemFile)
self.assertTrue(
np.all(cutout.getImage().array == cutoutFromBytes.getImage().array))
cutoutBytes = packageAlerts.streamCcdDataToBytes(cutoutCcdData)
with io.BytesIO(cutoutBytes) as bytesIO:
cutoutFromBytes = CCDData.read(bytesIO, format="fits")
self.assertFloatsAlmostEqual(cutoutCcdData.data, cutoutFromBytes.data)

def testMakeAlertDict(self):
"""Test stripping data from the various data products and into a
Expand All @@ -319,21 +367,24 @@ def testMakeAlertDict(self):
cutout = self.exposure.getCutout(sphPoint,
geom.Extent2I(self.cutoutSize,
self.cutoutSize))
cutputBytes = packageAlerts.makeCutoutBytes(cutout)
ccdCutout = packageAlerts.createCcdDataCutout(
cutout, sphPoint, cutout.getPhotoCalib())
cutoutBytes = packageAlerts.streamCcdDataToBytes(
ccdCutout)
objSources = self.diaSourceHistory.loc[srcIdx[0]]
alert = packageAlerts.makeAlertDict(
alertId,
diaSource,
self.diaObjects.loc[srcIdx[0]],
objSources,
cutout,
ccdCutout,
None)
self.assertEqual(len(alert), 9)

self.assertEqual(alert["alertId"], alertId)
self.assertEqual(alert["diaSource"], diaSource.to_dict())
self.assertEqual(alert["cutoutDifference"],
cutputBytes)
cutoutBytes)
self.assertEqual(alert["cutoutTemplate"],
None)

Expand Down Expand Up @@ -374,8 +425,10 @@ def testRun(self):
cutout = self.exposure.getCutout(sphPoint,
geom.Extent2I(self.cutoutSize,
self.cutoutSize))
ccdCutout = packageAlerts.createCcdDataCutout(
cutout, sphPoint, cutout.getPhotoCalib())
self.assertEqual(alert["cutoutDifference"],
packageAlerts.makeCutoutBytes(cutout))
packageAlerts.streamCcdDataToBytes(ccdCutout))

shutil.rmtree(tempdir)

Expand Down