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

DM-32835: Do injection source selection on sky instead of CCD #614

Open
wants to merge 4 commits into
base: main
Choose a base branch
from
Open

DM-32835: Do injection source selection on sky instead of CCD #614

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
0ed83e7
Don't call deprecated overload of psf.computeApertureFlux
jmeyers314 Dec 6, 2021
8f65f93
Do CCD overlap selection on the sky
jmeyers314 Dec 6, 2021
a4d5bd2
Elide addPixCoords when possible
jmeyers314 Dec 6, 2021
ad2107a
WCS area check no longer required
jmeyers314 Dec 6, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
59 changes: 33 additions & 26 deletions python/lsst/pipe/tasks/insertFakes.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@
import lsst.afw.math as afwMath
import lsst.pex.config as pexConfig
import lsst.pipe.base as pipeBase
import lsst.sphgeom as sphgeom

from lsst.pipe.base import CmdLineTask, PipelineTask, PipelineTaskConfig, PipelineTaskConnections
import lsst.pipe.base.connectionTypes as cT
Expand Down Expand Up @@ -71,8 +72,6 @@ def _add_fake_sources(exposure, objects, calibFluxRadius=12.0, logger=None):
fullBounds = galsim.BoundsI(bbox.minX, bbox.maxX, bbox.minY, bbox.maxY)
gsImg = galsim.Image(exposure.image.array, bounds=fullBounds)

pixScale = wcs.getPixelScale().asArcseconds()

for spt, gsObj in objects:
pt = wcs.skyToPixel(spt)
posd = galsim.PositionD(pt.x, pt.y)
Expand All @@ -83,15 +82,10 @@ def _add_fake_sources(exposure, objects, calibFluxRadius=12.0, logger=None):
mat = wcs.linearizePixelToSky(spt, geom.arcseconds).getMatrix()
gsWCS = galsim.JacobianWCS(mat[0, 0], mat[0, 1], mat[1, 0], mat[1, 1])

# This check is here because sometimes the WCS
# is multivalued and objects that should not be
# were being included.
gsPixScale = np.sqrt(gsWCS.pixelArea())
if gsPixScale < pixScale/2 or gsPixScale > pixScale*2:
continue

try:
psfArr = psf.computeKernelImage(pt).array
apCorr = psf.computeApertureFlux(calibFluxRadius, pt)
psfArr /= apCorr
except InvalidParameterError:
# Try mapping to nearest point contained in bbox.
contained_pt = Point2D(
Expand All @@ -108,6 +102,8 @@ def _add_fake_sources(exposure, objects, calibFluxRadius=12.0, logger=None):
# otherwise, try again with new point
try:
psfArr = psf.computeKernelImage(contained_pt).array
apCorr = psf.computeApertureFlux(calibFluxRadius, contained_pt)
psfArr /= apCorr
except InvalidParameterError:
if logger:
logger.infof(
Expand All @@ -116,8 +112,6 @@ def _add_fake_sources(exposure, objects, calibFluxRadius=12.0, logger=None):
)
continue

apCorr = psf.computeApertureFlux(calibFluxRadius)
psfArr /= apCorr
gsPSF = galsim.InterpolatedImage(galsim.Image(psfArr), wcs=gsWCS)

conv = galsim.Convolve(gsObj, gsPSF)
Expand Down Expand Up @@ -633,8 +627,6 @@ def run(self, fakeCat, image, wcs, photoCalib):

band = image.getFilterLabel().bandLabel
fakeCat = self._standardizeColumns(fakeCat, band)

fakeCat = self.addPixCoords(fakeCat, image)
fakeCat = self.trimFakeCat(fakeCat, image)

if len(fakeCat) > 0:
Expand Down Expand Up @@ -945,6 +937,9 @@ def processImagesForInsertion(self, fakeCat, wcs, psf, photoCalib, band, pixelSc

self.log.info("Processing %d fake images", len(fakeCat))

if 'x' not in fakeCat.columns:
fakeCat = self.addPixCoords(fakeCat, wcs=wcs)

for (imFile, sourceType, mag, x, y) in zip(fakeCat[band + "imFilename"].array,
fakeCat["sourceType"].array,
fakeCat['mag'].array,
Expand Down Expand Up @@ -991,22 +986,25 @@ def processImagesForInsertion(self, fakeCat, wcs, psf, photoCalib, band, pixelSc

return galImages, starImages

def addPixCoords(self, fakeCat, image):
def addPixCoords(self, fakeCat, image=None, wcs=None):

"""Add pixel coordinates to the catalog of fakes.

Parameters
----------
fakeCat : `pandas.core.frame.DataFrame`
The catalog of fake sources to be input
image : `lsst.afw.image.exposure.exposure.ExposureF`
image : `lsst.afw.image.exposure.exposure.ExposureF`, optional
The image into which the fake sources should be added
wcs : `lsst.afw.geom.SkyWcs`, optional
WCS to use to add fake sources

Returns
-------
fakeCat : `pandas.core.frame.DataFrame`
"""
wcs = image.getWcs()
if wcs is None:
wcs = image.getWcs()
ras = fakeCat['ra'].values
decs = fakeCat['dec'].values
xs, ys = wcs.skyToPixelArray(ras, decs)
Expand All @@ -1018,8 +1016,6 @@ def addPixCoords(self, fakeCat, image):
def trimFakeCat(self, fakeCat, image):
"""Trim the fake cat to about the size of the input image.

`fakeCat` must be processed with addPixCoords before using this method.

Parameters
----------
fakeCat : `pandas.core.frame.DataFrame`
Expand All @@ -1032,15 +1028,20 @@ def trimFakeCat(self, fakeCat, image):
fakeCat : `pandas.core.frame.DataFrame`
The original fakeCat trimmed to the area of the image
"""

bbox = Box2D(image.getBBox()).dilatedBy(self.config.trimBuffer)
xs = fakeCat["x"].values
ys = fakeCat["y"].values

isContained = xs >= bbox.minX
isContained &= xs <= bbox.maxX
isContained &= ys >= bbox.minY
isContained &= ys <= bbox.maxY
sphgeomCorners = []
for corner in bbox.getCorners():
spt = image.wcs.pixelToSky(corner)
sphgeomCorners.append(
sphgeom.UnitVector3d(
sphgeom.NormalizedAngle.fromRadians(spt[0].asRadians()),
sphgeom.Angle.fromRadians(spt[1].asRadians())
)
)
poly = sphgeom.ConvexPolygon(sphgeomCorners)
isContained = poly.contains(
fakeCat['ra'], fakeCat['dec']
)

return fakeCat[isContained]

Expand Down Expand Up @@ -1080,6 +1081,9 @@ def mkFakeGalsimGalaxies(self, fakeCat, band, photoCalib, pixelScale, psf, image

self.log.info("Making %d fake galaxy images", len(fakeCat))

if 'x' not in fakeCat.columns:
fakeCat = self.addPixCoords(fakeCat, image)

for (index, row) in fakeCat.iterrows():
xy = geom.Point2D(row["x"], row["y"])

Expand Down Expand Up @@ -1153,6 +1157,9 @@ def mkFakeStars(self, fakeCat, band, photoCalib, psf, image):

self.log.info("Making %d fake star images", len(fakeCat))

if 'x' not in fakeCat.columns:
fakeCat = self.addPixCoords(fakeCat, image)

for (index, row) in fakeCat.iterrows():
xy = geom.Point2D(row["x"], row["y"])

Expand Down