LSST Applications g0265f82a02+c6dfa2ddaf,g1162b98a3f+ffe7eabc7e,g2079a07aa2+1b2e822518,g2bbee38e9b+c6dfa2ddaf,g337abbeb29+c6dfa2ddaf,g36da64cc00+ea84795170,g3ddfee87b4+955a963fd8,g50ff169b8f+2eb0e556e8,g52b1c1532d+90ebb246c7,g555ede804d+955a963fd8,g591dd9f2cf+bac198a2cb,g5ec818987f+420292cfeb,g858d7b2824+d6c9a0a3b8,g876c692160+aabc49a3c3,g8a8a8dda67+90ebb246c7,g8cdfe0ae6a+4fd9e222a8,g99cad8db69+e6cd765486,g9ddcbc5298+a1346535a5,ga1e77700b3+df8f93165b,ga8c6da7877+acd47f83f4,gae46bcf261+c6dfa2ddaf,gb0e22166c9+8634eb87fb,gb3f2274832+12c8382528,gba4ed39666+1ac82b564f,gbb8dafda3b+0574160a1f,gbeb006f7da+dea2fbb49f,gc28159a63d+c6dfa2ddaf,gc86a011abf+d6c9a0a3b8,gcf0d15dbbd+955a963fd8,gdaeeff99f8+1cafcb7cd4,gdc0c513512+d6c9a0a3b8,ge79ae78c31+c6dfa2ddaf,geb67518f79+ba1859f325,gee10cc3b42+90ebb246c7,gf1cff7945b+d6c9a0a3b8,w.2024.13
LSST Data Management Base Package
Loading...
Searching...
No Matches
Classes | Functions | Variables
lsst.ip.diffim.utils Namespace Reference

Classes

class  DipoleTestImage
 

Functions

 showSourceSet (sSet, xy0=(0, 0), frame=0, ctype=afwDisplay.GREEN, symb="+", size=2)
 
 showKernelSpatialCells (maskedIm, kernelCellSet, showChi2=False, symb="o", ctype=None, ctypeUnused=None, ctypeBad=None, size=3, frame=None, title="Spatial Cells")
 
 showDiaSources (sources, exposure, isFlagged, isDipole, frame=None)
 
 showKernelCandidates (kernelCellSet, kernel, background, frame=None, showBadCandidates=True, resids=False, kernels=False)
 
 showKernelBasis (kernel, frame=None)
 
 plotKernelSpatialModel (kernel, kernelCellSet, showBadCandidates=True, numSample=128, keepPlots=True, maxCoeff=10)
 
 plotKernelCoefficients (spatialKernel, kernelCellSet, showBadCandidates=False, keepPlots=True)
 
 showKernelMosaic (bbox, kernel, nx=7, ny=None, frame=None, title=None, showCenter=True, showEllipticity=True)
 
 plotWhisker (results, newWcs)
 
 _sliceWidth (image, threshold, peaks, axis)
 
 getPsfFwhm (psf, average=True, position=None)
 
float evaluateMeanPsfFwhm (afwImage.Exposure exposure, float fwhmExposureBuffer, int fwhmExposureGrid)
 
afwImage.ImageD computeAveragePsf (afwImage.Exposure exposure, float psfExposureBuffer, int psfExposureGrid)
 
 detectTestSources (exposure, addMaskPlanes=None)
 
 makeFakeWcs ()
 
 makeTestImage (seed=5, nSrc=20, psfSize=2., noiseLevel=5., noiseSeed=6, fluxLevel=500., fluxRange=2., kernelSize=32, templateBorderSize=0, background=None, xSize=256, ySize=256, x0=12345, y0=67890, calibration=1., doApplyCalibration=False, xLoc=None, yLoc=None, flux=None, clearEdgeMask=False, addMaskPlanes=None)
 
 _makeTruthSchema ()
 
 _fillTruthCatalog (injectList)
 
 makeStats (badMaskPlanes=None)
 
 computeRobustStatistics (image, mask, statsCtrl, statistic=afwMath.MEANCLIP)
 
 computePSFNoiseEquivalentArea (psf)
 

Variables

bool keptPlots = False
 
 _LOG = getLogger(__name__)
 

Detailed Description

Support utilities for Measuring sources

Function Documentation

◆ _fillTruthCatalog()

lsst.ip.diffim.utils._fillTruthCatalog ( injectList)
protected
Add injected sources to the truth catalog.

Parameters
----------
injectList : `list` [`float`]
    Sources that were injected; tuples of (x, y, flux, size).

Returns
-------
catalog : `lsst.afw.table.SourceCatalog`
    Catalog with centroids and instFlux/instFluxErr values filled in and
    appropriate slots set.

Definition at line 1233 of file utils.py.

1233def _fillTruthCatalog(injectList):
1234 """Add injected sources to the truth catalog.
1235
1236 Parameters
1237 ----------
1238 injectList : `list` [`float`]
1239 Sources that were injected; tuples of (x, y, flux, size).
1240
1241 Returns
1242 -------
1243 catalog : `lsst.afw.table.SourceCatalog`
1244 Catalog with centroids and instFlux/instFluxErr values filled in and
1245 appropriate slots set.
1246 """
1247 keys, schema = _makeTruthSchema()
1248 catalog = afwTable.SourceCatalog(schema)
1249 catalog.reserve(len(injectList))
1250 for x, y, flux, size in injectList:
1251 record = catalog.addNew()
1252 keys["centroid"].set(record, geom.PointD(x, y))
1253 keys["instFlux"].set(record, flux)
1254 # Approximate injected errors
1255 keys["instFluxErr"].set(record, 20)
1256 # 5-sigma effective source width
1257 circle = afwGeom.Ellipse(afwGeom.ellipses.Axes(5*size, 5*size, 0), geom.Point2D(x, y))
1258 footprint = afwDetection.Footprint(afwGeom.SpanSet.fromShape(circle))
1259 footprint.addPeak(x, y, flux)
1260 record.setFootprint(footprint)
1261
1262 return catalog
1263
1264
An ellipse defined by an arbitrary BaseCore and a center point.
Definition Ellipse.h:51
daf::base::PropertySet * set
Definition fits.cc:931

◆ _makeTruthSchema()

lsst.ip.diffim.utils._makeTruthSchema ( )
protected
Make a schema for the truth catalog produced by `makeTestImage`.

Returns
-------
keys : `dict` [`str`]
    Fields added to the catalog, to make it easier to set them.
schema : `lsst.afw.table.Schema`
    Schema to use to make a "truth" SourceCatalog.
    Calib, Ap, and Psf flux slots all are set to ``truth_instFlux``.

Definition at line 1201 of file utils.py.

1201def _makeTruthSchema():
1202 """Make a schema for the truth catalog produced by `makeTestImage`.
1203
1204 Returns
1205 -------
1206 keys : `dict` [`str`]
1207 Fields added to the catalog, to make it easier to set them.
1208 schema : `lsst.afw.table.Schema`
1209 Schema to use to make a "truth" SourceCatalog.
1210 Calib, Ap, and Psf flux slots all are set to ``truth_instFlux``.
1211 """
1212 schema = afwTable.SourceTable.makeMinimalSchema()
1213 keys = {}
1214 # Don't use a FluxResultKey so we can manage the flux and err separately.
1215 keys["instFlux"] = schema.addField("truth_instFlux", type=np.float64,
1216 doc="true instFlux", units="count")
1217 keys["instFluxErr"] = schema.addField("truth_instFluxErr", type=np.float64,
1218 doc="true instFluxErr", units="count")
1219 keys["centroid"] = afwTable.Point2DKey.addFields(schema, "truth", "true simulated centroid", "pixel")
1220 schema.addField("truth_flag", "Flag", "truth flux failure flag.")
1221 # Add the flag fields a source selector would need.
1222 schema.addField("sky_source", "Flag", "testing flag.")
1223 schema.addField("base_PixelFlags_flag_interpolated", "Flag", "testing flag.")
1224 schema.addField("base_PixelFlags_flag_saturated", "Flag", "testing flag.")
1225 schema.addField("base_PixelFlags_flag_bad", "Flag", "testing flag.")
1226 schema.getAliasMap().set("slot_Centroid", "truth")
1227 schema.getAliasMap().set("slot_CalibFlux", "truth")
1228 schema.getAliasMap().set("slot_ApFlux", "truth")
1229 schema.getAliasMap().set("slot_PsfFlux", "truth")
1230 return keys, schema
1231
1232

◆ _sliceWidth()

lsst.ip.diffim.utils._sliceWidth ( image,
threshold,
peaks,
axis )
protected

Definition at line 856 of file utils.py.

856def _sliceWidth(image, threshold, peaks, axis):
857 vec = image.take(peaks[1 - axis], axis=axis)
858 low = np.interp(threshold, vec[:peaks[axis] + 1], np.arange(peaks[axis] + 1))
859 high = np.interp(threshold, vec[:peaks[axis] - 1:-1], np.arange(len(vec) - 1, peaks[axis] - 1, -1))
860 return high - low
861
862

◆ computeAveragePsf()

afwImage.ImageD lsst.ip.diffim.utils.computeAveragePsf ( afwImage.Exposure exposure,
float psfExposureBuffer,
int psfExposureGrid )
Get the average PSF by evaluating it on a grid within an exposure.

Parameters
----------
exposure : `~lsst.afw.image.Exposure`
    The exposure for which the average PSF is to be computed.
    The exposure must contain a `psf` attribute.
psfExposureBuffer : `float`
    Fractional buffer margin to be left out of all sides of the image
    during the construction of the grid to compute average PSF in an
    exposure.
psfExposureGrid : `int`
    Grid size to compute the average PSF in an exposure.

Returns
-------
psfImage : `~lsst.afw.image.Image`
    The average PSF across the exposure.

Raises
------
ValueError
    Raised if the PSF cannot be computed at any of the grid points.

See Also
--------
`evaluateMeanPsfFwhm`

Definition at line 957 of file utils.py.

958 psfExposureBuffer: float, psfExposureGrid: int) -> afwImage.ImageD:
959 """Get the average PSF by evaluating it on a grid within an exposure.
960
961 Parameters
962 ----------
963 exposure : `~lsst.afw.image.Exposure`
964 The exposure for which the average PSF is to be computed.
965 The exposure must contain a `psf` attribute.
966 psfExposureBuffer : `float`
967 Fractional buffer margin to be left out of all sides of the image
968 during the construction of the grid to compute average PSF in an
969 exposure.
970 psfExposureGrid : `int`
971 Grid size to compute the average PSF in an exposure.
972
973 Returns
974 -------
975 psfImage : `~lsst.afw.image.Image`
976 The average PSF across the exposure.
977
978 Raises
979 ------
980 ValueError
981 Raised if the PSF cannot be computed at any of the grid points.
982
983 See Also
984 --------
985 `evaluateMeanPsfFwhm`
986 """
987
988 psf = exposure.psf
989
990 bbox = exposure.getBBox()
991 xmax, ymax = bbox.getMax()
992 xmin, ymin = bbox.getMin()
993
994 xbuffer = psfExposureBuffer*(xmax-xmin)
995 ybuffer = psfExposureBuffer*(ymax-ymin)
996
997 nImg = 0
998 psfArray = None
999 for (x, y) in itertools.product(np.linspace(xmin+xbuffer, xmax-xbuffer, psfExposureGrid),
1000 np.linspace(ymin+ybuffer, ymax-ybuffer, psfExposureGrid)
1001 ):
1002 pos = geom.Point2D(x, y)
1003 try:
1004 singleImage = psf.computeKernelImage(pos)
1005 except InvalidParameterError:
1006 _LOG.debug("Unable to compute PSF image at position (%f, %f).", x, y)
1007 continue
1008
1009 if psfArray is None:
1010 psfArray = singleImage.array
1011 else:
1012 psfArray += singleImage.array
1013 nImg += 1
1014
1015 if psfArray is None:
1016 raise ValueError("Unable to compute PSF image at any position on the exposure.")
1017
1018 psfImage = afwImage.ImageD(psfArray/nImg)
1019 return psfImage
1020
1021

◆ computePSFNoiseEquivalentArea()

lsst.ip.diffim.utils.computePSFNoiseEquivalentArea ( psf)
Compute the noise equivalent area for an image psf

Parameters
----------
psf : `lsst.afw.detection.Psf`

Returns
-------
nea : `float`

Definition at line 1312 of file utils.py.

1312def computePSFNoiseEquivalentArea(psf):
1313 """Compute the noise equivalent area for an image psf
1314
1315 Parameters
1316 ----------
1317 psf : `lsst.afw.detection.Psf`
1318
1319 Returns
1320 -------
1321 nea : `float`
1322 """
1323 psfImg = psf.computeImage(psf.getAveragePosition())
1324 nea = 1./np.sum(psfImg.array**2)
1325 return nea

◆ computeRobustStatistics()

lsst.ip.diffim.utils.computeRobustStatistics ( image,
mask,
statsCtrl,
statistic = afwMath.MEANCLIP )
Calculate a robust mean of the variance plane of an exposure.

Parameters
----------
image : `lsst.afw.image.Image`
    Image or variance plane of an exposure to evaluate.
mask : `lsst.afw.image.Mask`
    Mask plane to use for excluding pixels.
statsCtrl : `lsst.afw.math.StatisticsControl`
    Statistics control object for configuring the calculation.
statistic : `lsst.afw.math.Property`, optional
    The type of statistic to compute. Typical values are
    ``afwMath.MEANCLIP`` or ``afwMath.STDEVCLIP``.

Returns
-------
value : `float`
    The result of the statistic calculated from the unflagged pixels.

Definition at line 1288 of file utils.py.

1288def computeRobustStatistics(image, mask, statsCtrl, statistic=afwMath.MEANCLIP):
1289 """Calculate a robust mean of the variance plane of an exposure.
1290
1291 Parameters
1292 ----------
1293 image : `lsst.afw.image.Image`
1294 Image or variance plane of an exposure to evaluate.
1295 mask : `lsst.afw.image.Mask`
1296 Mask plane to use for excluding pixels.
1297 statsCtrl : `lsst.afw.math.StatisticsControl`
1298 Statistics control object for configuring the calculation.
1299 statistic : `lsst.afw.math.Property`, optional
1300 The type of statistic to compute. Typical values are
1301 ``afwMath.MEANCLIP`` or ``afwMath.STDEVCLIP``.
1302
1303 Returns
1304 -------
1305 value : `float`
1306 The result of the statistic calculated from the unflagged pixels.
1307 """
1308 statObj = afwMath.makeStatistics(image, mask, statistic, statsCtrl)
1309 return statObj.getValue(statistic)
1310
1311
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
Definition Statistics.h:361

◆ detectTestSources()

lsst.ip.diffim.utils.detectTestSources ( exposure,
addMaskPlanes = None )
Minimal source detection wrapper suitable for unit tests.

Parameters
----------
exposure : `lsst.afw.image.Exposure`
    Exposure on which to run detection/measurement
    The exposure is modified in place to set the 'DETECTED' mask plane.
addMaskPlanes : `list` of `str`, optional
    Additional mask planes to add to the maskedImage of the exposure.

Returns
-------
selectSources
    Source catalog containing candidates

Definition at line 1022 of file utils.py.

1022def detectTestSources(exposure, addMaskPlanes=None):
1023 """Minimal source detection wrapper suitable for unit tests.
1024
1025 Parameters
1026 ----------
1027 exposure : `lsst.afw.image.Exposure`
1028 Exposure on which to run detection/measurement
1029 The exposure is modified in place to set the 'DETECTED' mask plane.
1030 addMaskPlanes : `list` of `str`, optional
1031 Additional mask planes to add to the maskedImage of the exposure.
1032
1033 Returns
1034 -------
1035 selectSources
1036 Source catalog containing candidates
1037 """
1038 if addMaskPlanes is None:
1039 # add empty streak mask plane in lieu of maskStreaksTask
1040 # And add empty INJECTED and INJECTED_TEMPLATE mask planes
1041 addMaskPlanes = ["STREAK", "INJECTED", "INJECTED_TEMPLATE"]
1042
1043 schema = afwTable.SourceTable.makeMinimalSchema()
1044 selectDetection = measAlg.SourceDetectionTask(schema=schema)
1045 selectMeasurement = measBase.SingleFrameMeasurementTask(schema=schema)
1046 table = afwTable.SourceTable.make(schema)
1047
1048 detRet = selectDetection.run(
1049 table=table,
1050 exposure=exposure,
1051 sigma=None, # The appropriate sigma is calculated from the PSF
1052 doSmooth=True
1053 )
1054 for mp in addMaskPlanes:
1055 exposure.mask.addMaskPlane(mp)
1056
1057 selectSources = detRet.sources
1058 selectMeasurement.run(measCat=selectSources, exposure=exposure)
1059
1060 return selectSources
1061
1062

◆ evaluateMeanPsfFwhm()

float lsst.ip.diffim.utils.evaluateMeanPsfFwhm ( afwImage.Exposure exposure,
float fwhmExposureBuffer,
int fwhmExposureGrid )
Get the mean PSF FWHM by evaluating it on a grid within an exposure.

Parameters
----------
exposure : `~lsst.afw.image.Exposure`
    The exposure for which the mean FWHM of the PSF is to be computed.
    The exposure must contain a `psf` attribute.
fwhmExposureBuffer : `float`
    Fractional buffer margin to be left out of all sides of the image
    during the construction of the grid to compute mean PSF FWHM in an
    exposure.
fwhmExposureGrid : `int`
    Grid size to compute the mean FWHM in an exposure.

Returns
-------
meanFwhm : `float`
    The mean PSF FWHM on the exposure.

Raises
------
ValueError
    Raised if the PSF cannot be computed at any of the grid points.

See Also
--------
`getPsfFwhm`
`computeAveragePsf`

Definition at line 897 of file utils.py.

898 fwhmExposureBuffer: float, fwhmExposureGrid: int) -> float:
899 """Get the mean PSF FWHM by evaluating it on a grid within an exposure.
900
901 Parameters
902 ----------
903 exposure : `~lsst.afw.image.Exposure`
904 The exposure for which the mean FWHM of the PSF is to be computed.
905 The exposure must contain a `psf` attribute.
906 fwhmExposureBuffer : `float`
907 Fractional buffer margin to be left out of all sides of the image
908 during the construction of the grid to compute mean PSF FWHM in an
909 exposure.
910 fwhmExposureGrid : `int`
911 Grid size to compute the mean FWHM in an exposure.
912
913 Returns
914 -------
915 meanFwhm : `float`
916 The mean PSF FWHM on the exposure.
917
918 Raises
919 ------
920 ValueError
921 Raised if the PSF cannot be computed at any of the grid points.
922
923 See Also
924 --------
925 `getPsfFwhm`
926 `computeAveragePsf`
927 """
928
929 psf = exposure.psf
930
931 bbox = exposure.getBBox()
932 xmax, ymax = bbox.getMax()
933 xmin, ymin = bbox.getMin()
934
935 xbuffer = fwhmExposureBuffer*(xmax-xmin)
936 ybuffer = fwhmExposureBuffer*(ymax-ymin)
937
938 width = []
939 for (x, y) in itertools.product(np.linspace(xmin+xbuffer, xmax-xbuffer, fwhmExposureGrid),
940 np.linspace(ymin+ybuffer, ymax-ybuffer, fwhmExposureGrid)
941 ):
942 pos = geom.Point2D(x, y)
943 try:
944 fwhm = getPsfFwhm(psf, average=True, position=pos)
945 except InvalidParameterError:
946 _LOG.debug("Unable to compute PSF FWHM at position (%f, %f).", x, y)
947 continue
948
949 width.append(fwhm)
950
951 if not width:
952 raise ValueError("Unable to compute PSF FWHM at any position on the exposure.")
953
954 return np.nanmean(width)
955
956

◆ getPsfFwhm()

lsst.ip.diffim.utils.getPsfFwhm ( psf,
average = True,
position = None )
Directly calculate the horizontal and vertical widths
of a PSF at half its maximum value.

Parameters
----------
psf : `~lsst.afw.detection.Psf`
    Point spread function (PSF) to evaluate.
average : `bool`, optional
    Set to return the average width over Y and X axes.
position : `~lsst.geom.Point2D`, optional
    The position at which to evaluate the PSF. If `None`, then the
    average position is used.

Returns
-------
psfSize : `float` | `tuple` [`float`]
    The FWHM of the PSF computed at its average position.
    Returns the widths along the Y and X axes,
    or the average of the two if `average` is set.

See Also
--------
evaluateMeanPsfFwhm

Definition at line 863 of file utils.py.

863def getPsfFwhm(psf, average=True, position=None):
864 """Directly calculate the horizontal and vertical widths
865 of a PSF at half its maximum value.
866
867 Parameters
868 ----------
869 psf : `~lsst.afw.detection.Psf`
870 Point spread function (PSF) to evaluate.
871 average : `bool`, optional
872 Set to return the average width over Y and X axes.
873 position : `~lsst.geom.Point2D`, optional
874 The position at which to evaluate the PSF. If `None`, then the
875 average position is used.
876
877 Returns
878 -------
879 psfSize : `float` | `tuple` [`float`]
880 The FWHM of the PSF computed at its average position.
881 Returns the widths along the Y and X axes,
882 or the average of the two if `average` is set.
883
884 See Also
885 --------
886 evaluateMeanPsfFwhm
887 """
888 if position is None:
889 position = psf.getAveragePosition()
890 image = psf.computeKernelImage(position).array
891 peak = psf.computePeak(position)
892 peakLocs = np.unravel_index(np.argmax(image), image.shape)
893 width = _sliceWidth(image, peak/2., peakLocs, axis=0), _sliceWidth(image, peak/2., peakLocs, axis=1)
894 return np.nanmean(width) if average else width
895
896

◆ makeFakeWcs()

lsst.ip.diffim.utils.makeFakeWcs ( )
Make a fake, affine Wcs.

Definition at line 1063 of file utils.py.

1063def makeFakeWcs():
1064 """Make a fake, affine Wcs.
1065 """
1066 crpix = geom.Point2D(123.45, 678.9)
1067 crval = geom.SpherePoint(0.1, 0.1, geom.degrees)
1068 cdMatrix = np.array([[5.19513851e-05, -2.81124812e-07],
1069 [-3.25186974e-07, -5.19112119e-05]])
1070 return afwGeom.makeSkyWcs(crpix, crval, cdMatrix)
1071
1072
Point in an unspecified spherical coordinate system.
Definition SpherePoint.h:57
std::shared_ptr< SkyWcs > makeSkyWcs(daf::base::PropertySet &metadata, bool strip=false)
Construct a SkyWcs from FITS keywords.
Definition SkyWcs.cc:521

◆ makeStats()

lsst.ip.diffim.utils.makeStats ( badMaskPlanes = None)
Create a statistics control for configuring calculations on images.

Parameters
----------
badMaskPlanes : `list` of `str`, optional
    List of mask planes to exclude from calculations.

Returns
-------
statsControl : ` lsst.afw.math.StatisticsControl`
    Statistics control object for configuring calculations on images.

Definition at line 1265 of file utils.py.

1265def makeStats(badMaskPlanes=None):
1266 """Create a statistics control for configuring calculations on images.
1267
1268 Parameters
1269 ----------
1270 badMaskPlanes : `list` of `str`, optional
1271 List of mask planes to exclude from calculations.
1272
1273 Returns
1274 -------
1275 statsControl : ` lsst.afw.math.StatisticsControl`
1276 Statistics control object for configuring calculations on images.
1277 """
1278 if badMaskPlanes is None:
1279 badMaskPlanes = ("INTRP", "EDGE", "DETECTED", "SAT", "CR",
1280 "BAD", "NO_DATA", "DETECTED_NEGATIVE")
1281 statsControl = afwMath.StatisticsControl()
1282 statsControl.setNumSigmaClip(3.)
1283 statsControl.setNumIter(3)
1284 statsControl.setAndMask(afwImage.Mask.getPlaneBitMask(badMaskPlanes))
1285 return statsControl
1286
1287
Pass parameters to a Statistics object.
Definition Statistics.h:83

◆ makeTestImage()

lsst.ip.diffim.utils.makeTestImage ( seed = 5,
nSrc = 20,
psfSize = 2.,
noiseLevel = 5.,
noiseSeed = 6,
fluxLevel = 500.,
fluxRange = 2.,
kernelSize = 32,
templateBorderSize = 0,
background = None,
xSize = 256,
ySize = 256,
x0 = 12345,
y0 = 67890,
calibration = 1.,
doApplyCalibration = False,
xLoc = None,
yLoc = None,
flux = None,
clearEdgeMask = False,
addMaskPlanes = None )
Make a reproduceable PSF-convolved exposure for testing.

Parameters
----------
seed : `int`, optional
    Seed value to initialize the random number generator for sources.
nSrc : `int`, optional
    Number of sources to simulate.
psfSize : `float`, optional
    Width of the PSF of the simulated sources, in pixels.
noiseLevel : `float`, optional
    Standard deviation of the noise to add to each pixel.
noiseSeed : `int`, optional
    Seed value to initialize the random number generator for noise.
fluxLevel : `float`, optional
    Reference flux of the simulated sources.
fluxRange : `float`, optional
    Range in flux amplitude of the simulated sources.
kernelSize : `int`, optional
    Size in pixels of the kernel for simulating sources.
templateBorderSize : `int`, optional
    Size in pixels of the image border used to pad the image.
background : `lsst.afw.math.Chebyshev1Function2D`, optional
    Optional background to add to the output image.
xSize, ySize : `int`, optional
    Size in pixels of the simulated image.
x0, y0 : `int`, optional
    Origin of the image.
calibration : `float`, optional
    Conversion factor between instFlux and nJy.
doApplyCalibration : `bool`, optional
    Apply the photometric calibration and return the image in nJy?
xLoc, yLoc : `list` of `float`, optional
    User-specified coordinates of the simulated sources.
    If specified, must have length equal to ``nSrc``
flux : `list` of `float`, optional
    User-specified fluxes of the simulated sources.
    If specified, must have length equal to ``nSrc``
clearEdgeMask : `bool`, optional
    Clear the "EDGE" mask plane after source detection.
addMaskPlanes : `list` of `str`, optional
    Mask plane names to add to the image.

Returns
-------
modelExposure : `lsst.afw.image.Exposure`
    The model image, with the mask and variance planes. The DETECTED
    plane is filled in for the injected source footprints.
sourceCat : `lsst.afw.table.SourceCatalog`
    Catalog of sources inserted in the model image.

Raises
------
ValueError
    If `xloc`, `yloc`, or `flux` are supplied with inconsistant lengths.

Definition at line 1073 of file utils.py.

1088 ):
1089 """Make a reproduceable PSF-convolved exposure for testing.
1090
1091 Parameters
1092 ----------
1093 seed : `int`, optional
1094 Seed value to initialize the random number generator for sources.
1095 nSrc : `int`, optional
1096 Number of sources to simulate.
1097 psfSize : `float`, optional
1098 Width of the PSF of the simulated sources, in pixels.
1099 noiseLevel : `float`, optional
1100 Standard deviation of the noise to add to each pixel.
1101 noiseSeed : `int`, optional
1102 Seed value to initialize the random number generator for noise.
1103 fluxLevel : `float`, optional
1104 Reference flux of the simulated sources.
1105 fluxRange : `float`, optional
1106 Range in flux amplitude of the simulated sources.
1107 kernelSize : `int`, optional
1108 Size in pixels of the kernel for simulating sources.
1109 templateBorderSize : `int`, optional
1110 Size in pixels of the image border used to pad the image.
1111 background : `lsst.afw.math.Chebyshev1Function2D`, optional
1112 Optional background to add to the output image.
1113 xSize, ySize : `int`, optional
1114 Size in pixels of the simulated image.
1115 x0, y0 : `int`, optional
1116 Origin of the image.
1117 calibration : `float`, optional
1118 Conversion factor between instFlux and nJy.
1119 doApplyCalibration : `bool`, optional
1120 Apply the photometric calibration and return the image in nJy?
1121 xLoc, yLoc : `list` of `float`, optional
1122 User-specified coordinates of the simulated sources.
1123 If specified, must have length equal to ``nSrc``
1124 flux : `list` of `float`, optional
1125 User-specified fluxes of the simulated sources.
1126 If specified, must have length equal to ``nSrc``
1127 clearEdgeMask : `bool`, optional
1128 Clear the "EDGE" mask plane after source detection.
1129 addMaskPlanes : `list` of `str`, optional
1130 Mask plane names to add to the image.
1131
1132 Returns
1133 -------
1134 modelExposure : `lsst.afw.image.Exposure`
1135 The model image, with the mask and variance planes. The DETECTED
1136 plane is filled in for the injected source footprints.
1137 sourceCat : `lsst.afw.table.SourceCatalog`
1138 Catalog of sources inserted in the model image.
1139
1140 Raises
1141 ------
1142 ValueError
1143 If `xloc`, `yloc`, or `flux` are supplied with inconsistant lengths.
1144 """
1145 # Distance from the inner edge of the bounding box to avoid placing test
1146 # sources in the model images.
1147 bufferSize = kernelSize/2 + templateBorderSize + 1
1148
1149 bbox = geom.Box2I(geom.Point2I(x0, y0), geom.Extent2I(xSize, ySize))
1150 if templateBorderSize > 0:
1151 bbox.grow(templateBorderSize)
1152
1153 rng = np.random.RandomState(seed)
1154 rngNoise = np.random.RandomState(noiseSeed)
1155 x0, y0 = bbox.getBegin()
1156 xSize, ySize = bbox.getDimensions()
1157 if xLoc is None:
1158 xLoc = rng.rand(nSrc)*(xSize - 2*bufferSize) + bufferSize + x0
1159 else:
1160 if len(xLoc) != nSrc:
1161 raise ValueError("xLoc must have length equal to nSrc. %f supplied vs %f", len(xLoc), nSrc)
1162 if yLoc is None:
1163 yLoc = rng.rand(nSrc)*(ySize - 2*bufferSize) + bufferSize + y0
1164 else:
1165 if len(yLoc) != nSrc:
1166 raise ValueError("yLoc must have length equal to nSrc. %f supplied vs %f", len(yLoc), nSrc)
1167
1168 if flux is None:
1169 flux = (rng.rand(nSrc)*(fluxRange - 1.) + 1.)*fluxLevel
1170 else:
1171 if len(flux) != nSrc:
1172 raise ValueError("flux must have length equal to nSrc. %f supplied vs %f", len(flux), nSrc)
1173 sigmas = [psfSize for src in range(nSrc)]
1174 injectList = list(zip(xLoc, yLoc, flux, sigmas))
1175 skyLevel = 0
1176 # Don't use the built in poisson noise: it modifies the global state of numpy random
1177 modelExposure = plantSources(bbox, kernelSize, skyLevel, injectList, addPoissonNoise=False)
1178 modelExposure.setWcs(makeFakeWcs())
1179 noise = rngNoise.randn(ySize, xSize)*noiseLevel
1180 noise -= np.mean(noise)
1181 modelExposure.variance.array = np.sqrt(np.abs(modelExposure.image.array)) + noiseLevel**2
1182 modelExposure.image.array += noise
1183
1184 # Run source detection to set up the mask plane
1185 detectTestSources(modelExposure, addMaskPlanes=addMaskPlanes)
1186 if clearEdgeMask:
1187 modelExposure.mask &= ~modelExposure.mask.getPlaneBitMask("EDGE")
1188 modelExposure.setPhotoCalib(afwImage.PhotoCalib(calibration, 0., bbox))
1189 if background is not None:
1190 modelExposure.image += background
1191 modelExposure.maskedImage /= calibration
1192 modelExposure.info.setId(seed)
1193 if doApplyCalibration:
1194 modelExposure.maskedImage = modelExposure.photoCalib.calibrateImage(modelExposure.maskedImage)
1195
1196 truth = _fillTruthCatalog(injectList)
1197
1198 return modelExposure, truth
1199
1200
The photometric calibration of an exposure.
Definition PhotoCalib.h:114
An integer coordinate rectangle.
Definition Box.h:55

◆ plotKernelCoefficients()

lsst.ip.diffim.utils.plotKernelCoefficients ( spatialKernel,
kernelCellSet,
showBadCandidates = False,
keepPlots = True )
Plot the individual kernel candidate and the spatial kernel solution coefficients.

Parameters
----------

spatialKernel : `lsst.afw.math.LinearCombinationKernel`
    The spatial spatialKernel solution model which is a spatially varying linear combination
    of the spatialKernel basis functions.
    Typically returned by `lsst.ip.diffim.SpatialKernelSolution.getSolutionPair()`.

kernelCellSet : `lsst.afw.math.SpatialCellSet`
    The spatial cells that was used for solution for the spatialKernel. They contain the
    local solutions of the AL kernel for the selected sources.

showBadCandidates : `bool`, optional
    If True, plot the coefficient values for kernel candidates where the solution was marked
    bad by the numerical algorithm. Defaults to False.

keepPlots: `bool`, optional
    If True, sets ``plt.show()`` to be called before the task terminates, so that the plots
    can be explored interactively. Defaults to True.

Notes
-----
This function produces 3 figures per image subtraction operation.
* A grid plot of the local solutions. Each grid cell corresponds to a proportional area in
  the image. In each cell, local kernel solution coefficients are plotted of kernel candidates (color)
  that fall into this area as a function of the kernel basis function number.
* A grid plot of the spatial solution. Each grid cell corresponds to a proportional area in
  the image. In each cell, the spatial solution coefficients are evaluated for the center of the cell.
* Histogram of the local solution coefficients. Red line marks the spatial solution value at
  center of the image.

This function is called if ``lsst.ip.diffim.psfMatch.plotKernelCoefficients==True`` in lsstDebug. This
function was implemented as part of DM-17825.

Definition at line 432 of file utils.py.

432def plotKernelCoefficients(spatialKernel, kernelCellSet, showBadCandidates=False, keepPlots=True):
433 """Plot the individual kernel candidate and the spatial kernel solution coefficients.
434
435 Parameters
436 ----------
437
438 spatialKernel : `lsst.afw.math.LinearCombinationKernel`
439 The spatial spatialKernel solution model which is a spatially varying linear combination
440 of the spatialKernel basis functions.
441 Typically returned by `lsst.ip.diffim.SpatialKernelSolution.getSolutionPair()`.
442
443 kernelCellSet : `lsst.afw.math.SpatialCellSet`
444 The spatial cells that was used for solution for the spatialKernel. They contain the
445 local solutions of the AL kernel for the selected sources.
446
447 showBadCandidates : `bool`, optional
448 If True, plot the coefficient values for kernel candidates where the solution was marked
449 bad by the numerical algorithm. Defaults to False.
450
451 keepPlots: `bool`, optional
452 If True, sets ``plt.show()`` to be called before the task terminates, so that the plots
453 can be explored interactively. Defaults to True.
454
455 Notes
456 -----
457 This function produces 3 figures per image subtraction operation.
458 * A grid plot of the local solutions. Each grid cell corresponds to a proportional area in
459 the image. In each cell, local kernel solution coefficients are plotted of kernel candidates (color)
460 that fall into this area as a function of the kernel basis function number.
461 * A grid plot of the spatial solution. Each grid cell corresponds to a proportional area in
462 the image. In each cell, the spatial solution coefficients are evaluated for the center of the cell.
463 * Histogram of the local solution coefficients. Red line marks the spatial solution value at
464 center of the image.
465
466 This function is called if ``lsst.ip.diffim.psfMatch.plotKernelCoefficients==True`` in lsstDebug. This
467 function was implemented as part of DM-17825.
468 """
469 try:
470 import matplotlib.pyplot as plt
471 except ImportError as e:
472 print("Unable to import matplotlib: %s" % e)
473 return
474
475 # Image dimensions
476 imgBBox = kernelCellSet.getBBox()
477 x0 = imgBBox.getBeginX()
478 y0 = imgBBox.getBeginY()
479 wImage = imgBBox.getWidth()
480 hImage = imgBBox.getHeight()
481 imgCenterX = imgBBox.getCenterX()
482 imgCenterY = imgBBox.getCenterY()
483
484 # Plot the local solutions
485 # ----
486
487 # Grid size
488 nX = 8
489 nY = 8
490 wCell = wImage / nX
491 hCell = hImage / nY
492
493 fig = plt.figure()
494 fig.suptitle("Kernel candidate parameters on an image grid")
495 arrAx = fig.subplots(nrows=nY, ncols=nX, sharex=True, sharey=True, gridspec_kw=dict(
496 wspace=0, hspace=0))
497
498 # Bottom left panel is for bottom left part of the image
499 arrAx = arrAx[::-1, :]
500
501 allParams = []
502 for cell in kernelCellSet.getCellList():
503 cellBBox = geom.Box2D(cell.getBBox())
504 # Determine which panel this spatial cell belongs to
505 iX = int((cellBBox.getCenterX() - x0)//wCell)
506 iY = int((cellBBox.getCenterY() - y0)//hCell)
507
508 for cand in cell.begin(False):
509 try:
510 kernel = cand.getKernel(cand.ORIG)
511 except Exception:
512 continue
513
514 if not showBadCandidates and cand.isBad():
515 continue
516
517 nKernelParams = kernel.getNKernelParameters()
518 kernelParams = np.array(kernel.getKernelParameters())
519 allParams.append(kernelParams)
520
521 if cand.isBad():
522 color = 'red'
523 else:
524 color = None
525 arrAx[iY, iX].plot(np.arange(nKernelParams), kernelParams, '.-',
526 color=color, drawstyle='steps-mid', linewidth=0.1)
527 for ax in arrAx.ravel():
528 ax.grid(True, axis='y')
529
530 # Plot histogram of the local parameters and the global solution at the image center
531 # ----
532
533 spatialFuncs = spatialKernel.getSpatialFunctionList()
534 nKernelParams = spatialKernel.getNKernelParameters()
535 nX = 8
536 fig = plt.figure()
537 fig.suptitle("Hist. of parameters marked with spatial solution at img center")
538 arrAx = fig.subplots(nrows=int(nKernelParams//nX)+1, ncols=nX)
539 arrAx = arrAx[::-1, :]
540 allParams = np.array(allParams)
541 for k in range(nKernelParams):
542 ax = arrAx.ravel()[k]
543 ax.hist(allParams[:, k], bins=20, edgecolor='black')
544 ax.set_xlabel('P{}'.format(k))
545 valueParam = spatialFuncs[k](imgCenterX, imgCenterY)
546 ax.axvline(x=valueParam, color='red')
547 ax.text(0.1, 0.9, '{:.1f}'.format(valueParam),
548 transform=ax.transAxes, backgroundcolor='lightsteelblue')
549
550 # Plot grid of the spatial solution
551 # ----
552
553 nX = 8
554 nY = 8
555 wCell = wImage / nX
556 hCell = hImage / nY
557 x0 += wCell / 2
558 y0 += hCell / 2
559
560 fig = plt.figure()
561 fig.suptitle("Spatial solution of kernel parameters on an image grid")
562 arrAx = fig.subplots(nrows=nY, ncols=nX, sharex=True, sharey=True, gridspec_kw=dict(
563 wspace=0, hspace=0))
564 arrAx = arrAx[::-1, :]
565 kernelParams = np.zeros(nKernelParams, dtype=float)
566
567 for iX in range(nX):
568 for iY in range(nY):
569 x = x0 + iX * wCell
570 y = y0 + iY * hCell
571 # Evaluate the spatial solution functions for this x,y location
572 kernelParams = [f(x, y) for f in spatialFuncs]
573 arrAx[iY, iX].plot(np.arange(nKernelParams), kernelParams, '.-', drawstyle='steps-mid')
574 arrAx[iY, iX].grid(True, axis='y')
575
576 global keptPlots
577 if keepPlots and not keptPlots:
578 # Keep plots open when done
579 def show():
580 print("%s: Please close plots when done." % __name__)
581 try:
582 plt.show()
583 except Exception:
584 pass
585 print("Plots closed, exiting...")
586 import atexit
587 atexit.register(show)
588 keptPlots = True
589
590
A floating-point coordinate rectangle geometry.
Definition Box.h:413

◆ plotKernelSpatialModel()

lsst.ip.diffim.utils.plotKernelSpatialModel ( kernel,
kernelCellSet,
showBadCandidates = True,
numSample = 128,
keepPlots = True,
maxCoeff = 10 )
Plot the Kernel spatial model.

Definition at line 282 of file utils.py.

283 numSample=128, keepPlots=True, maxCoeff=10):
284 """Plot the Kernel spatial model.
285 """
286 try:
287 import matplotlib.pyplot as plt
288 import matplotlib.colors
289 except ImportError as e:
290 print("Unable to import numpy and matplotlib: %s" % e)
291 return
292
293 x0 = kernelCellSet.getBBox().getBeginX()
294 y0 = kernelCellSet.getBBox().getBeginY()
295
296 candPos = list()
297 candFits = list()
298 badPos = list()
299 badFits = list()
300 candAmps = list()
301 badAmps = list()
302 for cell in kernelCellSet.getCellList():
303 for cand in cell.begin(False):
304 if not showBadCandidates and cand.isBad():
305 continue
306 candCenter = geom.PointD(cand.getXCenter(), cand.getYCenter())
307 try:
308 im = cand.getTemplateMaskedImage()
309 except Exception:
310 continue
311
312 targetFits = badFits if cand.isBad() else candFits
313 targetPos = badPos if cand.isBad() else candPos
314 targetAmps = badAmps if cand.isBad() else candAmps
315
316 # compare original and spatial kernel coefficients
317 kp0 = np.array(cand.getKernel(diffimLib.KernelCandidateF.ORIG).getKernelParameters())
318 amp = cand.getCandidateRating()
319
320 targetFits = badFits if cand.isBad() else candFits
321 targetPos = badPos if cand.isBad() else candPos
322 targetAmps = badAmps if cand.isBad() else candAmps
323
324 targetFits.append(kp0)
325 targetPos.append(candCenter)
326 targetAmps.append(amp)
327
328 xGood = np.array([pos.getX() for pos in candPos]) - x0
329 yGood = np.array([pos.getY() for pos in candPos]) - y0
330 zGood = np.array(candFits)
331
332 xBad = np.array([pos.getX() for pos in badPos]) - x0
333 yBad = np.array([pos.getY() for pos in badPos]) - y0
334 zBad = np.array(badFits)
335 numBad = len(badPos)
336
337 xRange = np.linspace(0, kernelCellSet.getBBox().getWidth(), num=numSample)
338 yRange = np.linspace(0, kernelCellSet.getBBox().getHeight(), num=numSample)
339
340 if maxCoeff:
341 maxCoeff = min(maxCoeff, kernel.getNKernelParameters())
342 else:
343 maxCoeff = kernel.getNKernelParameters()
344
345 for k in range(maxCoeff):
346 func = kernel.getSpatialFunction(k)
347 dfGood = zGood[:, k] - np.array([func(pos.getX(), pos.getY()) for pos in candPos])
348 yMin = dfGood.min()
349 yMax = dfGood.max()
350 if numBad > 0:
351 dfBad = zBad[:, k] - np.array([func(pos.getX(), pos.getY()) for pos in badPos])
352 # Can really screw up the range...
353 yMin = min([yMin, dfBad.min()])
354 yMax = max([yMax, dfBad.max()])
355 yMin -= 0.05*(yMax - yMin)
356 yMax += 0.05*(yMax - yMin)
357
358 fRange = np.ndarray((len(xRange), len(yRange)))
359 for j, yVal in enumerate(yRange):
360 for i, xVal in enumerate(xRange):
361 fRange[j][i] = func(xVal, yVal)
362
363 fig = plt.figure(k)
364
365 fig.clf()
366 try:
367 fig.canvas._tkcanvas._root().lift() # == Tk's raise, but raise is a python reserved word
368 except Exception: # protect against API changes
369 pass
370
371 fig.suptitle('Kernel component %d' % k)
372
373 # LL
374 ax = fig.add_axes((0.1, 0.05, 0.35, 0.35))
375 vmin = fRange.min() # - 0.05*np.fabs(fRange.min())
376 vmax = fRange.max() # + 0.05*np.fabs(fRange.max())
377 norm = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
378 im = ax.imshow(fRange, aspect='auto', norm=norm,
379 extent=[0, kernelCellSet.getBBox().getWidth() - 1,
380 0, kernelCellSet.getBBox().getHeight() - 1])
381 ax.set_title('Spatial polynomial')
382 plt.colorbar(im, orientation='horizontal', ticks=[vmin, vmax])
383
384 # UL
385 ax = fig.add_axes((0.1, 0.55, 0.35, 0.35))
386 ax.plot(-2.5*np.log10(candAmps), zGood[:, k], 'b+')
387 if numBad > 0:
388 ax.plot(-2.5*np.log10(badAmps), zBad[:, k], 'r+')
389 ax.set_title("Basis Coefficients")
390 ax.set_xlabel("Instr mag")
391 ax.set_ylabel("Coeff")
392
393 # LR
394 ax = fig.add_axes((0.55, 0.05, 0.35, 0.35))
395 ax.set_autoscale_on(False)
396 ax.set_xbound(lower=0, upper=kernelCellSet.getBBox().getHeight())
397 ax.set_ybound(lower=yMin, upper=yMax)
398 ax.plot(yGood, dfGood, 'b+')
399 if numBad > 0:
400 ax.plot(yBad, dfBad, 'r+')
401 ax.axhline(0.0)
402 ax.set_title('dCoeff (indiv-spatial) vs. y')
403
404 # UR
405 ax = fig.add_axes((0.55, 0.55, 0.35, 0.35))
406 ax.set_autoscale_on(False)
407 ax.set_xbound(lower=0, upper=kernelCellSet.getBBox().getWidth())
408 ax.set_ybound(lower=yMin, upper=yMax)
409 ax.plot(xGood, dfGood, 'b+')
410 if numBad > 0:
411 ax.plot(xBad, dfBad, 'r+')
412 ax.axhline(0.0)
413 ax.set_title('dCoeff (indiv-spatial) vs. x')
414
415 fig.show()
416
417 global keptPlots
418 if keepPlots and not keptPlots:
419 # Keep plots open when done
420 def show():
421 print("%s: Please close plots when done." % __name__)
422 try:
423 plt.show()
424 except Exception:
425 pass
426 print("Plots closed, exiting...")
427 import atexit
428 atexit.register(show)
429 keptPlots = True
430
431
int min
int max

◆ plotWhisker()

lsst.ip.diffim.utils.plotWhisker ( results,
newWcs )
Plot whisker diagram of astromeric offsets between results.matches.

Definition at line 674 of file utils.py.

674def plotWhisker(results, newWcs):
675 """Plot whisker diagram of astromeric offsets between results.matches.
676 """
677 refCoordKey = results.matches[0].first.getTable().getCoordKey()
678 inCentroidKey = results.matches[0].second.getTable().getCentroidSlot().getMeasKey()
679 positions = [m.first.get(refCoordKey) for m in results.matches]
680 residuals = [m.first.get(refCoordKey).getOffsetFrom(
681 newWcs.pixelToSky(m.second.get(inCentroidKey))) for
682 m in results.matches]
683 import matplotlib.pyplot as plt
684 fig = plt.figure()
685 sp = fig.add_subplot(1, 1, 0)
686 xpos = [x[0].asDegrees() for x in positions]
687 ypos = [x[1].asDegrees() for x in positions]
688 xpos.append(0.02*(max(xpos) - min(xpos)) + min(xpos))
689 ypos.append(0.98*(max(ypos) - min(ypos)) + min(ypos))
690 xidxs = np.isfinite(xpos)
691 yidxs = np.isfinite(ypos)
692 X = np.asarray(xpos)[xidxs]
693 Y = np.asarray(ypos)[yidxs]
694 distance = [x[1].asArcseconds() for x in residuals]
695 distance.append(0.2)
696 distance = np.asarray(distance)[xidxs]
697 # NOTE: This assumes that the bearing is measured positive from +RA through North.
698 # From the documentation this is not clear.
699 bearing = [x[0].asRadians() for x in residuals]
700 bearing.append(0)
701 bearing = np.asarray(bearing)[xidxs]
702 U = (distance*np.cos(bearing))
703 V = (distance*np.sin(bearing))
704 sp.quiver(X, Y, U, V)
705 sp.set_title("WCS Residual")
706 plt.show()
707
708

◆ showDiaSources()

lsst.ip.diffim.utils.showDiaSources ( sources,
exposure,
isFlagged,
isDipole,
frame = None )
Display Dia Sources.

Definition at line 110 of file utils.py.

110def showDiaSources(sources, exposure, isFlagged, isDipole, frame=None):
111 """Display Dia Sources.
112 """
113 #
114 # Show us the ccandidates
115 #
116 # Too many mask planes in diffims
117 disp = afwDisplay.Display(frame=frame)
118 for plane in ("BAD", "CR", "EDGE", "INTERPOlATED", "INTRP", "SAT", "SATURATED"):
119 disp.setMaskPlaneColor(plane, color="ignore")
120
121 mos = afwDisplay.utils.Mosaic()
122 for i in range(len(sources)):
123 source = sources[i]
124 badFlag = isFlagged[i]
125 dipoleFlag = isDipole[i]
126 bbox = source.getFootprint().getBBox()
127 stamp = exposure.Factory(exposure, bbox, True)
128 im = afwDisplay.utils.Mosaic(gutter=1, background=0, mode="x")
129 im.append(stamp.getMaskedImage())
130 lab = "%.1f,%.1f:" % (source.getX(), source.getY())
131 if badFlag:
132 ctype = afwDisplay.RED
133 lab += "BAD"
134 if dipoleFlag:
135 ctype = afwDisplay.YELLOW
136 lab += "DIPOLE"
137 if not badFlag and not dipoleFlag:
138 ctype = afwDisplay.GREEN
139 lab += "OK"
140 mos.append(im.makeMosaic(), lab, ctype)
141 title = "Dia Sources"
142 mosaicImage = mos.makeMosaic(display=disp, title=title)
143 return mosaicImage
144
145

◆ showKernelBasis()

lsst.ip.diffim.utils.showKernelBasis ( kernel,
frame = None )
Display a Kernel's basis images.

Definition at line 264 of file utils.py.

264def showKernelBasis(kernel, frame=None):
265 """Display a Kernel's basis images.
266 """
267 mos = afwDisplay.utils.Mosaic()
268
269 for k in kernel.getKernelList():
270 im = afwImage.ImageD(k.getDimensions())
271 k.computeImage(im, False)
272 mos.append(im)
273
274 disp = afwDisplay.Display(frame=frame)
275 mos.makeMosaic(display=disp, title="Kernel Basis Images")
276
277 return mos
278

◆ showKernelCandidates()

lsst.ip.diffim.utils.showKernelCandidates ( kernelCellSet,
kernel,
background,
frame = None,
showBadCandidates = True,
resids = False,
kernels = False )
Display the Kernel candidates.

If kernel is provided include spatial model and residuals;
If chi is True, generate a plot of residuals/sqrt(variance), i.e. chi.

Definition at line 146 of file utils.py.

147 resids=False, kernels=False):
148 """Display the Kernel candidates.
149
150 If kernel is provided include spatial model and residuals;
151 If chi is True, generate a plot of residuals/sqrt(variance), i.e. chi.
152 """
153 #
154 # Show us the ccandidates
155 #
156 if kernels:
157 mos = afwDisplay.utils.Mosaic(gutter=5, background=0)
158 else:
159 mos = afwDisplay.utils.Mosaic(gutter=5, background=-1)
160 #
161 candidateCenters = []
162 candidateCentersBad = []
163 candidateIndex = 0
164 for cell in kernelCellSet.getCellList():
165 for cand in cell.begin(False): # include bad candidates
166 # Original difference image; if does not exist, skip candidate
167 try:
168 resid = cand.getDifferenceImage(diffimLib.KernelCandidateF.ORIG)
169 except Exception:
170 continue
171
172 rchi2 = cand.getChi2()
173 if rchi2 > 1e100:
174 rchi2 = np.nan
175
176 if not showBadCandidates and cand.isBad():
177 continue
178
179 im_resid = afwDisplay.utils.Mosaic(gutter=1, background=-0.5, mode="x")
180
181 try:
182 im = cand.getScienceMaskedImage()
183 im = im.Factory(im, True)
184 im.setXY0(cand.getScienceMaskedImage().getXY0())
185 except Exception:
186 continue
187 if (not resids and not kernels):
188 im_resid.append(im.Factory(im, True))
189 try:
190 im = cand.getTemplateMaskedImage()
191 im = im.Factory(im, True)
192 im.setXY0(cand.getTemplateMaskedImage().getXY0())
193 except Exception:
194 continue
195 if (not resids and not kernels):
196 im_resid.append(im.Factory(im, True))
197
198 # Difference image with original basis
199 if resids:
200 var = resid.variance
201 var = var.Factory(var, True)
202 np.sqrt(var.array, var.array) # inplace sqrt
203 resid = resid.image
204 resid /= var
205 bbox = kernel.shrinkBBox(resid.getBBox())
206 resid = resid.Factory(resid, bbox, deep=True)
207 elif kernels:
208 kim = cand.getKernelImage(diffimLib.KernelCandidateF.ORIG).convertF()
209 resid = kim.Factory(kim, True)
210 im_resid.append(resid)
211
212 # residuals using spatial model
213 ski = afwImage.ImageD(kernel.getDimensions())
214 kernel.computeImage(ski, False, int(cand.getXCenter()), int(cand.getYCenter()))
215 sk = afwMath.FixedKernel(ski)
216 sbg = 0.0
217 if background:
218 sbg = background(int(cand.getXCenter()), int(cand.getYCenter()))
219 sresid = cand.getDifferenceImage(sk, sbg)
220 resid = sresid
221 if resids:
222 resid = sresid.image
223 resid /= var
224 bbox = kernel.shrinkBBox(resid.getBBox())
225 resid = resid.Factory(resid, bbox, deep=True)
226 elif kernels:
227 kim = ski.convertF()
228 resid = kim.Factory(kim, True)
229 im_resid.append(resid)
230
231 im = im_resid.makeMosaic()
232
233 lab = "%d chi^2 %.1f" % (cand.getId(), rchi2)
234 ctype = afwDisplay.RED if cand.isBad() else afwDisplay.GREEN
235
236 mos.append(im, lab, ctype)
237
238 if False and np.isnan(rchi2):
239 disp = afwDisplay.Display(frame=1)
240 disp.mtv(cand.getScienceMaskedImage.image, title="candidate")
241 print("rating", cand.getCandidateRating())
242
243 im = cand.getScienceMaskedImage()
244 center = (candidateIndex, cand.getXCenter() - im.getX0(), cand.getYCenter() - im.getY0())
245 candidateIndex += 1
246 if cand.isBad():
247 candidateCentersBad.append(center)
248 else:
249 candidateCenters.append(center)
250
251 if resids:
252 title = "chi Diffim"
253 elif kernels:
254 title = "Kernels"
255 else:
256 title = "Candidates & residuals"
257
258 disp = afwDisplay.Display(frame=frame)
259 mosaicImage = mos.makeMosaic(display=disp, title=title)
260
261 return mosaicImage
262
263
A kernel created from an Image.
Definition Kernel.h:471

◆ showKernelMosaic()

lsst.ip.diffim.utils.showKernelMosaic ( bbox,
kernel,
nx = 7,
ny = None,
frame = None,
title = None,
showCenter = True,
showEllipticity = True )
Show a mosaic of Kernel images.

Definition at line 591 of file utils.py.

592 showCenter=True, showEllipticity=True):
593 """Show a mosaic of Kernel images.
594 """
595 mos = afwDisplay.utils.Mosaic()
596
597 x0 = bbox.getBeginX()
598 y0 = bbox.getBeginY()
599 width = bbox.getWidth()
600 height = bbox.getHeight()
601
602 if not ny:
603 ny = int(nx*float(height)/width + 0.5)
604 if not ny:
605 ny = 1
606
607 schema = afwTable.SourceTable.makeMinimalSchema()
608 centroidName = "base_SdssCentroid"
609 shapeName = "base_SdssShape"
610 control = measBase.SdssCentroidControl()
611 schema.getAliasMap().set("slot_Centroid", centroidName)
612 schema.getAliasMap().set("slot_Centroid_flag", centroidName + "_flag")
613 centroider = measBase.SdssCentroidAlgorithm(control, centroidName, schema)
614 sdssShape = measBase.SdssShapeControl()
615 shaper = measBase.SdssShapeAlgorithm(sdssShape, shapeName, schema)
616 table = afwTable.SourceTable.make(schema)
617 table.defineCentroid(centroidName)
618 table.defineShape(shapeName)
619
620 centers = []
621 shapes = []
622 for iy in range(ny):
623 for ix in range(nx):
624 x = int(ix*(width - 1)/(nx - 1)) + x0
625 y = int(iy*(height - 1)/(ny - 1)) + y0
626
627 im = afwImage.ImageD(kernel.getDimensions())
628 ksum = kernel.computeImage(im, False, x, y)
629 lab = "Kernel(%d,%d)=%.2f" % (x, y, ksum) if False else ""
630 mos.append(im, lab)
631
632 # SdssCentroidAlgorithm.measure requires an exposure of floats
634
635 w, h = im.getWidth(), im.getHeight()
636 centerX = im.getX0() + w//2
637 centerY = im.getY0() + h//2
638 src = table.makeRecord()
639 spans = afwGeom.SpanSet(exp.getBBox())
640 foot = afwDet.Footprint(spans)
641 foot.addPeak(centerX, centerY, 1)
642 src.setFootprint(foot)
643
644 try: # The centroider requires a psf, so this will fail if none is attached to exp
645 centroider.measure(src, exp)
646 centers.append((src.getX(), src.getY()))
647
648 shaper.measure(src, exp)
649 shapes.append((src.getIxx(), src.getIxy(), src.getIyy()))
650 except Exception:
651 pass
652
653 disp = afwDisplay.Display(frame=frame)
654 mos.makeMosaic(display=disp, title=title if title else "Model Kernel", mode=nx)
655
656 if centers and frame is not None:
657 disp = afwDisplay.Display(frame=frame)
658 i = 0
659 with disp.Buffering():
660 for cen, shape in zip(centers, shapes):
661 bbox = mos.getBBox(i)
662 i += 1
663 xc, yc = cen[0] + bbox.getMinX(), cen[1] + bbox.getMinY()
664 if showCenter:
665 disp.dot("+", xc, yc, ctype=afwDisplay.BLUE)
666
667 if showEllipticity:
668 ixx, ixy, iyy = shape
669 disp.dot("@:%g,%g,%g" % (ixx, ixy, iyy), xc, yc, ctype=afwDisplay.RED)
670
671 return mos
672
673
Class to describe the properties of a detected object from an image.
Definition Footprint.h:63
A compact representation of a collection of pixels.
Definition SpanSet.h:78
MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > * makeMaskedImage(typename std::shared_ptr< Image< ImagePixelT > > image, typename std::shared_ptr< Mask< MaskPixelT > > mask=Mask< MaskPixelT >(), typename std::shared_ptr< Image< VariancePixelT > > variance=Image< VariancePixelT >())
A function to return a MaskedImage of the correct type (cf.
std::shared_ptr< Exposure< ImagePixelT, MaskPixelT, VariancePixelT > > makeExposure(MaskedImage< ImagePixelT, MaskPixelT, VariancePixelT > &mimage, std::shared_ptr< geom::SkyWcs const > wcs=std::shared_ptr< geom::SkyWcs const >())
A function to return an Exposure of the correct type (cf.
Definition Exposure.h:484

◆ showKernelSpatialCells()

lsst.ip.diffim.utils.showKernelSpatialCells ( maskedIm,
kernelCellSet,
showChi2 = False,
symb = "o",
ctype = None,
ctypeUnused = None,
ctypeBad = None,
size = 3,
frame = None,
title = "Spatial Cells" )
Show the SpatialCells.

If symb is something that display.dot understands (e.g. "o"), the top
nMaxPerCell candidates will be indicated with that symbol, using ctype
and size.

Definition at line 71 of file utils.py.

73 frame=None, title="Spatial Cells"):
74 """Show the SpatialCells.
75
76 If symb is something that display.dot understands (e.g. "o"), the top
77 nMaxPerCell candidates will be indicated with that symbol, using ctype
78 and size.
79 """
80 disp = afwDisplay.Display(frame=frame)
81 disp.mtv(maskedIm, title=title)
82 with disp.Buffering():
83 origin = [-maskedIm.getX0(), -maskedIm.getY0()]
84 for cell in kernelCellSet.getCellList():
85 afwDisplay.utils.drawBBox(cell.getBBox(), origin=origin, display=disp)
86
87 goodies = ctypeBad is None
88 for cand in cell.begin(goodies):
89 xc, yc = cand.getXCenter() + origin[0], cand.getYCenter() + origin[1]
90 if cand.getStatus() == afwMath.SpatialCellCandidate.BAD:
91 color = ctypeBad
92 elif cand.getStatus() == afwMath.SpatialCellCandidate.GOOD:
93 color = ctype
94 elif cand.getStatus() == afwMath.SpatialCellCandidate.UNKNOWN:
95 color = ctypeUnused
96 else:
97 continue
98
99 if color:
100 disp.dot(symb, xc, yc, ctype=color, size=size)
101
102 if showChi2:
103 rchi2 = cand.getChi2()
104 if rchi2 > 1e100:
105 rchi2 = np.nan
106 disp.dot("%d %.1f" % (cand.getId(), rchi2),
107 xc - size, yc - size - 4, ctype=color, size=size)
108
109

◆ showSourceSet()

lsst.ip.diffim.utils.showSourceSet ( sSet,
xy0 = (0, 0),
frame = 0,
ctype = afwDisplay.GREEN,
symb = "+",
size = 2 )
Draw the (XAstrom, YAstrom) positions of a set of Sources.

Image has the given XY0.

Definition at line 51 of file utils.py.

51def showSourceSet(sSet, xy0=(0, 0), frame=0, ctype=afwDisplay.GREEN, symb="+", size=2):
52 """Draw the (XAstrom, YAstrom) positions of a set of Sources.
53
54 Image has the given XY0.
55 """
56 disp = afwDisplay.afwDisplay(frame=frame)
57 with disp.Buffering():
58 for s in sSet:
59 xc, yc = s.getXAstrom() - xy0[0], s.getYAstrom() - xy0[1]
60
61 if symb == "id":
62 disp.dot(str(s.getId()), xc, yc, ctype=ctype, size=size)
63 else:
64 disp.dot(symb, xc, yc, ctype=ctype, size=size)
65
66
67# Kernel display utilities
68#
69
70

Variable Documentation

◆ _LOG

lsst.ip.diffim.utils._LOG = getLogger(__name__)
protected

Definition at line 48 of file utils.py.

◆ keptPlots

bool lsst.ip.diffim.utils.keptPlots = False

Definition at line 46 of file utils.py.