Skip to content

Commit

Permalink
Merge pull request #185 from lsst/tickets/DM-41956
Browse files Browse the repository at this point in the history 
DM-41956: Exclude diaObjects not contained in the image bounding box.
  • Loading branch information
@isullivan
isullivan authored Dec 9, 2023
2 parents 4934932 + 4f67843 commit e8a2b74
Show file tree
Hide file tree
Showing 4 changed files with 171 additions and 72 deletions.
52 changes: 50 additions & 2 deletions python/lsst/ap/association/diaPipe.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@
"DiaPipelineTask",
"DiaPipelineConnections")

import numpy as np
import pandas as pd

from lsst.daf.base import DateTime
Expand Down Expand Up @@ -264,6 +265,13 @@ class DiaPipelineConfig(pipeBase.PipelineTaskConfig,
doc="Write out associated DiaSources, DiaForcedSources, and DiaObjects, "
"formatted following the Science Data Model.",
)
imagePixelMargin = pexConfig.RangeField(
dtype=int,
default=10,
min=0,
doc="Pad the image by this many pixels before removing off-image "
"diaObjects for association.",
)
idGenerator = DetectorVisitIdGeneratorConfig.make_field()

def setDefaults(self):
Expand Down Expand Up @@ -374,9 +382,14 @@ def run(self,
"""
# Load the DiaObjects and DiaSource history.
loaderResult = self.diaCatalogLoader.run(diffIm, self.apdb)
if len(loaderResult.diaObjects) > 0:
diaObjects = self.purgeDiaObjects(diffIm.getBBox(), diffIm.getWcs(), loaderResult.diaObjects,
buffer=self.config.imagePixelMargin)
else:
diaObjects = loaderResult.diaObjects

# Associate new DiaSources with existing DiaObjects.
assocResults = self.associator.run(diaSourceTable, loaderResult.diaObjects,
assocResults = self.associator.run(diaSourceTable, diaObjects,
exposure_time=diffIm.visitInfo.exposureTime)

if self.config.doSolarSystemAssociation:
Expand Down Expand Up @@ -419,7 +432,7 @@ def run(self,

# Append new DiaObjects and DiaSources to their previous history.
diaObjects = pd.concat(
[loaderResult.diaObjects,
[diaObjects,
createResults.newDiaObjects.set_index("diaObjectId", drop=False)],
sort=True)
if self.testDataFrameIndex(diaObjects):
Expand Down Expand Up @@ -634,3 +647,38 @@ def _add_association_meta_data(self,
self.metadata.add('numNewDiaObjects', nNewDiaObjects)
self.metadata.add('numTotalSolarSystemObjects', nTotalSsObjects)
self.metadata.add('numAssociatedSsObjects', nAssociatedSsObjects)

def purgeDiaObjects(self, bbox, wcs, diaObjCat, buffer=0):
"""Drop diaObjects that are outside the exposure bounding box.
Parameters
----------
bbox : `lsst.geom.Box2I`
Bounding box of the exposure.
wcs : `lsst.afw.geom.SkyWcs`
Coordinate system definition (wcs) for the exposure.
diaObjCat : `pandas.DataFrame`
DiaObjects loaded from the Apdb.
buffer : `int`, optional
Width, in pixels, to pad the exposure bounding box.
Returns
-------
diaObjCat : `pandas.DataFrame`
DiaObjects loaded from the Apdb, restricted to the exposure
bounding box.
"""
try:
bbox.grow(buffer)
raVals = diaObjCat.ra.to_numpy()
decVals = diaObjCat.dec.to_numpy()
xVals, yVals = wcs.skyToPixelArray(raVals, decVals, degrees=True)
selector = bbox.contains(xVals, yVals)
nPurged = np.sum(~selector)
if nPurged > 0:
diaObjCat = diaObjCat[selector].copy()
self.log.info(f"Dropped {nPurged} diaObjects that were outside the bbox "
f"leaving {len(diaObjCat)} in the catalog")
except Exception as e:
self.log.warning("Error attempting to check diaObject history: %s", e)
return diaObjCat
52 changes: 52 additions & 0 deletions tests/test_diaPipe.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@
import lsst.afw.image as afwImage
import lsst.afw.table as afwTable
from lsst.pipe.base.testUtils import assertValidOutput
from utils_tests import makeExposure, makeDiaObjects
import lsst.utils.tests
import lsst.utils.timer
from unittest.mock import patch, Mock, MagicMock, DEFAULT
Expand Down Expand Up @@ -181,6 +182,57 @@ def test_createDiaObjects(self):
result.newDiaObjects["diaObjectId"].to_numpy(),
result.diaSources["diaSourceId"].to_numpy())))

def test_purgeDiaObjects(self):
"""Remove diaOjects that are outside an image's bounding box.
"""

config = self._makeDefaultConfig(doPackageAlerts=False)
task = DiaPipelineTask(config=config)
exposure = makeExposure(False, False)
nObj0 = 20

# Create diaObjects
diaObjects = makeDiaObjects(nObj0, exposure)
# Shrink the bounding box so that some of the diaObjects will be outside
bbox = exposure.getBBox()
size = np.minimum(bbox.getHeight(), bbox.getWidth())
bbox.grow(-size//4)
exposureCut = exposure[bbox]
sizeCut = np.minimum(bbox.getHeight(), bbox.getWidth())
buffer = 10
bbox.grow(buffer)

def check_diaObjects(bbox, wcs, diaObjects):
raVals = diaObjects.ra.to_numpy()
decVals = diaObjects.dec.to_numpy()
xVals, yVals = wcs.skyToPixelArray(raVals, decVals, degrees=True)
selector = bbox.contains(xVals, yVals)
return selector

selector0 = check_diaObjects(bbox, exposureCut.getWcs(), diaObjects)
nIn0 = np.count_nonzero(selector0)
nOut0 = np.count_nonzero(~selector0)
self.assertEqual(nObj0, nIn0 + nOut0)

diaObjects1 = task.purgeDiaObjects(exposureCut.getBBox(), exposureCut.getWcs(), diaObjects,
buffer=buffer)
# Verify that the bounding box was not changed
sizeCheck = np.minimum(exposureCut.getBBox().getHeight(), exposureCut.getBBox().getWidth())
self.assertEqual(sizeCut, sizeCheck)
selector1 = check_diaObjects(bbox, exposureCut.getWcs(), diaObjects1)
nIn1 = np.count_nonzero(selector1)
nOut1 = np.count_nonzero(~selector1)
nObj1 = len(diaObjects1)
self.assertEqual(nObj1, nIn0)
# Verify that not all diaObjects were removed
self.assertGreater(nObj1, 0)
# Check that some diaObjects were removed
self.assertLess(nObj1, nObj0)
# Verify that no objects outside the bounding box remain
self.assertEqual(nOut1, 0)
# Verify that no objects inside the bounding box were removed
self.assertEqual(nIn1, nIn0)


class MemoryTester(lsst.utils.tests.MemoryTestCase):
pass
Expand Down
74 changes: 4 additions & 70 deletions tests/test_loadDiaCatalogs.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,15 +25,11 @@
import unittest
import yaml

from lsst.afw.cameraGeom.testUtils import DetectorWrapper
import lsst.afw.geom as afwGeom
import lsst.afw.image as afwImage
from lsst.ap.association import LoadDiaCatalogsTask, LoadDiaCatalogsConfig
import lsst.daf.base as dafBase
from lsst.dax.apdb import ApdbSql, ApdbSqlConfig, ApdbTables
from lsst.utils import getPackageDir
import lsst.utils.tests
import utils_tests
from utils_tests import makeExposure, makeDiaObjects, makeDiaSources, makeDiaForcedSources


def _data_file_name(basename, module_name):
Expand All @@ -55,68 +51,6 @@ def _data_file_name(basename, module_name):
return os.path.join(getPackageDir(module_name), "data", basename)


def makeExposure(flipX=False, flipY=False):
"""Create an exposure and flip the x or y (or both) coordinates.
Returns bounding boxes that are right or left handed around the bounding
polygon.
Parameters
----------
flipX : `bool`
Flip the x coordinate in the WCS.
flipY : `bool`
Flip the y coordinate in the WCS.
Returns
-------
exposure : `lsst.afw.image.Exposure`
Exposure with a valid bounding box and wcs.
"""
metadata = dafBase.PropertySet()

metadata.set("SIMPLE", "T")
metadata.set("BITPIX", -32)
metadata.set("NAXIS", 2)
metadata.set("NAXIS1", 1024)
metadata.set("NAXIS2", 1153)
metadata.set("RADECSYS", 'FK5')
metadata.set("EQUINOX", 2000.)

metadata.setDouble("CRVAL1", 215.604025685476)
metadata.setDouble("CRVAL2", 53.1595451514076)
metadata.setDouble("CRPIX1", 1109.99981456774)
metadata.setDouble("CRPIX2", 560.018167811613)
metadata.set("CTYPE1", 'RA---SIN')
metadata.set("CTYPE2", 'DEC--SIN')

xFlip = 1
if flipX:
xFlip = -1
yFlip = 1
if flipY:
yFlip = -1
metadata.setDouble("CD1_1", xFlip * 5.10808596133527E-05)
metadata.setDouble("CD1_2", yFlip * 1.85579539217196E-07)
metadata.setDouble("CD2_2", yFlip * -5.10281493481982E-05)
metadata.setDouble("CD2_1", xFlip * -8.27440751733828E-07)

wcs = afwGeom.makeSkyWcs(metadata)
exposure = afwImage.makeExposure(
afwImage.makeMaskedImageFromArrays(np.ones((1024, 1153))), wcs)
detector = DetectorWrapper(id=23, bbox=exposure.getBBox()).detector
visit = afwImage.VisitInfo(
exposureTime=200.,
date=dafBase.DateTime("2014-05-13T17:00:00.000000000",
dafBase.DateTime.Timescale.TAI))
exposure.info.id = 1234
exposure.setDetector(detector)
exposure.info.setVisitInfo(visit)
exposure.setFilter(afwImage.FilterLabel(band='g'))

return exposure


class TestLoadDiaCatalogs(unittest.TestCase):

def setUp(self):
Expand All @@ -135,12 +69,12 @@ def setUp(self):

self.exposure = makeExposure(False, False)

self.diaObjects = utils_tests.makeDiaObjects(20, self.exposure)
self.diaSources = utils_tests.makeDiaSources(
self.diaObjects = makeDiaObjects(20, self.exposure)
self.diaSources = makeDiaSources(
100,
self.diaObjects["diaObjectId"].to_numpy(),
self.exposure)
self.diaForcedSources = utils_tests.makeDiaForcedSources(
self.diaForcedSources = makeDiaForcedSources(
200,
self.diaObjects["diaObjectId"].to_numpy(),
self.exposure)
Expand Down
65 changes: 65 additions & 0 deletions tests/utils_tests.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,9 @@
import pandas as pd
import numpy as np

from lsst.afw.cameraGeom.testUtils import DetectorWrapper
import lsst.afw.geom as afwGeom
import lsst.afw.image as afwImage
import lsst.daf.base as dafBase
import lsst.geom

Expand Down Expand Up @@ -165,3 +168,65 @@ def makeDiaForcedSources(nForcedSources, diaObjectIds, exposure, randomizeObject
"flags": 0})

return pd.DataFrame(data=data)


def makeExposure(flipX=False, flipY=False):
"""Create an exposure and flip the x or y (or both) coordinates.
Returns bounding boxes that are right or left handed around the bounding
polygon.
Parameters
----------
flipX : `bool`
Flip the x coordinate in the WCS.
flipY : `bool`
Flip the y coordinate in the WCS.
Returns
-------
exposure : `lsst.afw.image.Exposure`
Exposure with a valid bounding box and wcs.
"""
metadata = dafBase.PropertySet()

metadata.set("SIMPLE", "T")
metadata.set("BITPIX", -32)
metadata.set("NAXIS", 2)
metadata.set("NAXIS1", 1024)
metadata.set("NAXIS2", 1153)
metadata.set("RADECSYS", 'FK5')
metadata.set("EQUINOX", 2000.)

metadata.setDouble("CRVAL1", 215.604025685476)
metadata.setDouble("CRVAL2", 53.1595451514076)
metadata.setDouble("CRPIX1", 1109.99981456774)
metadata.setDouble("CRPIX2", 560.018167811613)
metadata.set("CTYPE1", 'RA---SIN')
metadata.set("CTYPE2", 'DEC--SIN')

xFlip = 1
if flipX:
xFlip = -1
yFlip = 1
if flipY:
yFlip = -1
metadata.setDouble("CD1_1", xFlip * 5.10808596133527E-05)
metadata.setDouble("CD1_2", yFlip * 1.85579539217196E-07)
metadata.setDouble("CD2_2", yFlip * -5.10281493481982E-05)
metadata.setDouble("CD2_1", xFlip * -8.27440751733828E-07)

wcs = afwGeom.makeSkyWcs(metadata)
exposure = afwImage.makeExposure(
afwImage.makeMaskedImageFromArrays(np.ones((1024, 1153))), wcs)
detector = DetectorWrapper(id=23, bbox=exposure.getBBox()).detector
visit = afwImage.VisitInfo(
exposureTime=200.,
date=dafBase.DateTime("2014-05-13T17:00:00.000000000",
dafBase.DateTime.Timescale.TAI))
exposure.info.id = 1234
exposure.setDetector(detector)
exposure.info.setVisitInfo(visit)
exposure.setFilter(afwImage.FilterLabel(band='g'))

return exposure

0 comments on commit e8a2b74

Please sign in to comment.