LSST Applications  21.0.0-147-g0e635eb1+1acddb5be5,22.0.0+052faf71bd,22.0.0+1ea9a8b2b2,22.0.0+6312710a6c,22.0.0+729191ecac,22.0.0+7589c3a021,22.0.0+9f079a9461,22.0.1-1-g7d6de66+b8044ec9de,22.0.1-1-g87000a6+536b1ee016,22.0.1-1-g8e32f31+6312710a6c,22.0.1-10-gd060f87+016f7cdc03,22.0.1-12-g9c3108e+df145f6f68,22.0.1-16-g314fa6d+c825727ab8,22.0.1-19-g93a5c75+d23f2fb6d8,22.0.1-19-gb93eaa13+aab3ef7709,22.0.1-2-g8ef0a89+b8044ec9de,22.0.1-2-g92698f7+9f079a9461,22.0.1-2-ga9b0f51+052faf71bd,22.0.1-2-gac51dbf+052faf71bd,22.0.1-2-gb66926d+6312710a6c,22.0.1-2-gcb770ba+09e3807989,22.0.1-20-g32debb5+b8044ec9de,22.0.1-23-gc2439a9a+fb0756638e,22.0.1-3-g496fd5d+09117f784f,22.0.1-3-g59f966b+1e6ba2c031,22.0.1-3-g849a1b8+f8b568069f,22.0.1-3-gaaec9c0+c5c846a8b1,22.0.1-32-g5ddfab5d3+60ce4897b0,22.0.1-4-g037fbe1+64e601228d,22.0.1-4-g8623105+b8044ec9de,22.0.1-5-g096abc9+d18c45d440,22.0.1-5-g15c806e+57f5c03693,22.0.1-7-gba73697+57f5c03693,master-g6e05de7fdc+c1283a92b8,master-g72cdda8301+729191ecac,w.2021.39
LSST Data Management Base Package
cModelDisplay.py
Go to the documentation of this file.
1 __all__ = ("displayReconstructedCmodel", "displayReconstrucedCmodelMpl",
2  "buildCModelImages", "reconstructCModel")
3 
4 import numpy as np
5 import matplotlib.pyplot as plt
6 
7 import lsst.afw.image as afwImage
8 from lsst.geom import Point2D
9 import lsst.meas.modelfit as measMod
10 import lsst.shapelet as shapelet
11 
12 
13 def displayReconstructedCmodel(exposure, record, config, display="mpl"):
14  """ Display an image, the Cmodel, and residuals
15 
16  Parameters
17  ----------
18  exposure : `lsst.afw.image.Exposure`
19  Exposure object that contains the source which was modeled
20  record : `lsst.afw.table.SourceRecord`
21  Record object which contains the measurements made on the source
22  config : `lsst.meas.modelfit.CModel(SingleFrame/Forced)Config`
23  Configuration object of the CModel plugin used in the measurement
24  process
25  display : `str`, optional
26  Display in which to render the data, may be mpl for matplotlib or afw
27  for afwDisplay, defaults to mpl
28 
29  Raises
30  ------
31  `NotImplementedError`
32  If the display backend specified is not implemented
33  """
34 
35  if display == "mpl":
36  displayReconstrucedCmodelMpl(exposure, record, config)
37  elif display == "afw":
38  raise NotImplementedError("The afw display backend is not yet "
39  "implemented")
40  else:
41  raise NotImplementedError("The display backend '{}' is not a supported"
42  "backend".format(display))
43 
44 
45 def displayReconstrucedCmodelMpl(exposure, record, config):
46  """ Display an image, the Cmodel, and residuals using matplotlib
47 
48  Parameters
49  ----------
50  exposure : `lsst.afw.image.Exposure`
51  Exposure object that contains the source which was modeled
52  record : `lsst.afw.table.SourceRecord`
53  Record object which contains the measurements made on the source
54  config : `lsst.meas.modelfit.CModel(SingleFrame/Forced)Config`
55  Configuration object of the CModel plugin used in the measurement
56  process
57  """
58 
59  subImage, devIm, expIm, jointIm = buildCModelImages(exposure, record,
60  config)
61  # Get the min and max values for the sub image to use for scaling
62  subImMin = subImage.array.min()
63  subImMax = subImage.array.max()
64 
65  # Calculate the scaling to use for the residual images
66  devResidual = subImage.array-devIm.array
67  expResidual = subImage.array-expIm.array
68  jointResidual = subImage.array-jointIm.array
69  residualList = (devResidual, expResidual, jointResidual)
70 
71  differences = [(x.max()-x.max()) for x in residualList]
72  maxDifferenceIndex = np.argmax(differences)
73  resImMin = residualList[maxDifferenceIndex].min()
74  resImMax = residualList[maxDifferenceIndex].max()
75 
76  # Build the image figure
77  fig, axes = plt.subplots(3, 3, sharex='col', sharey='row', figsize=(8, 5))
78  fig.subplots_adjust(left=0.25, right=0.8)
79  lCBar = fig.add_axes([0.1, 0.15, 0.05, 0.7])
80  rCBar = fig.add_axes([0.85, 0.15, 0.05, 0.7])
81 
82  # Populate just the exposures in the appropriate places
83  for i in range(3):
84  axes[i, 0].imshow(subImage.array, vmin=subImMin, vmax=subImMax)
85 
86  # Populate dev panels
87  axes[0, 1].imshow(devIm.array, vmin=subImMin, vmax=subImMax)
88  axes[0, 2].imshow(devResidual, vmin=resImMin, vmax=resImMax, cmap="BrBG")
89 
90  # Populate exp panels
91  axes[1, 1].imshow(expIm.array, vmin=subImMin, vmax=subImMax)
92  axes[1, 2].imshow(expResidual, vmin=resImMin, vmax=resImMax, cmap="BrBG")
93 
94  # Populate joint panels
95  axes[2, 1].imshow(jointIm.array, vmin=subImMin, vmax=subImMax)
96  axes[2, 2].imshow(jointResidual, vmin=resImMin, vmax=resImMax, cmap="BrBG")
97 
98  axes[0, 0].set_title("Image")
99  axes[0, 1].set_title("Model")
100  axes[0, 2].set_title('Residuals')
101 
102  axes[0, 0].set_ylabel("Dev")
103  axes[1, 0].set_ylabel("Exp")
104  axes[2, 0].set_ylabel("Joint")
105 
106  fig.colorbar(axes[0, 0].get_images()[0], lCBar)
107  lCBar.yaxis.set_ticks_position('left')
108  fig.colorbar(axes[maxDifferenceIndex, 2].get_images()[0], rCBar)
109 
110  plt.show()
111 
112 
113 def buildCModelImages(exposure, record, config):
114  """ Create Images out of the CModel for the given record
115 
116  Parameters
117  ----------
118  exposure : `lsst.afw.image.Exposure`
119  Exposure object that contains the source which was modeled
120  record : `lsst.afw.table.SourceRecord`
121  Record object which contains the measurements made on the source
122  config : `lsst.meas.modelfit.CModel(SingleFrame/Forced)Config`
123  Configuration object of the CModel plugin used in the measurement
124  process
125 
126  Returns
127  -------
128  subImage : `lsst.afw.image.ImageD`
129  Sub image of the original data taken from a region defined by the
130  bounding box of the footprint for the object defined in the given
131  source record
132  devIm : `lsst.afw.image.ImageD`
133  Image created from the dev component of the CModel for the supplied
134  record at the same pixel locations as subImage
135  expIm: `lsst.afw.image.ImageD`
136  Image created from the exp component of the CModel for the supplied
137  record at the same pixel locations as subImage
138  jointIm :
139  Image created from the joint fit of the dev and exp components of the
140  CModel for the supplied record at the same pixel locations as subImage
141  """
142 
143  dev, exp, jointDev, jointExp = reconstructCModel(exposure, record, config)
144  # Get exposure cutout
145  footBBox = record.getFootprint().getBBox()
146  subImage = afwImage.ImageD(exposure.getImage().convertD(), footBBox)
147 
148  # Build the psf
149  shapeletPsfKey = shapelet.MultiShapeletFunctionKey(
150  record.schema[config.psfName])
151  psfApprox = record.get(shapeletPsfKey)
152 
153  # Build the dev Image from the shapelet function
154  devIm = afwImage.ImageD(footBBox)
155  dev = dev.convolve(psfApprox)
156  dev.evaluate().addToImage(devIm)
157 
158  # Build the exp image from the shapelet function
159  expIm = afwImage.ImageD(footBBox)
160  exp = exp.convolve(psfApprox)
161  exp.evaluate().addToImage(expIm)
162 
163  # Build the joint image from the shapelet function
164  jointIm = afwImage.ImageD(footBBox)
165  jointDev = jointDev.convolve(psfApprox)
166  jointExp = jointExp.convolve(psfApprox)
167  jointDev.evaluate().addToImage(jointIm)
168  jointExp.evaluate().addToImage(jointIm)
169 
170  return subImage, devIm, expIm, jointIm
171 
172 
173 def reconstructCModel(exposure, record, config):
174  """ Reconstruct the CModel for the given record
175 
176  Parameters
177  ----------
178  exposure : `lsst.afw.image.Exposure`
179  Exposure object that contains the source which was modeled
180  record : `lsst.afw.table.SourceRecord`
181  Record object which contains the measurements made on the source
182  config : `lsst.meas.modelfit.CModel(SingleFrame/Forced)Config`
183  Configuration object of the CModel plugin used in the measurement
184  process
185 
186  Returns
187  -------
188  devShapelet : `lsst.shapelet.MultiShapeletFunction`
189  Multi-component shapelet model of the dev component of CModel
190  expShapelet : `lsst.shapelet.MultiShapeletFunction`
191  Multi-component shapelet model fo the exp component of CModel
192  devJointShapelet : `lsst.shapelet.MultiShapeletFunction`
193  Multi-component shapelet model of the dev component of CModel jointly
194  fit with the exp component
195  expJointShapelet : `lsst.shapelet.MultiShapeletFunction`
196  Multi-component shapelet model of the exp component of Cmodel jointly
197  fit with the dev component
198  """
199 
200  # build a unit system transformation object
201  center = record.getCentroid()
202  position = exposure.getWcs().pixelToSky(center)
203  measSys = measMod.UnitSystem(exposure)
204  approxFlux = record.get("base_PsfFlux_instFlux")
205  fitSys = measMod.UnitSystem(position, exposure.getPhotoCalib(), approxFlux)
206  fitSysToMeasSys = measMod.LocalUnitTransform(Point2D(0, 0), fitSys, measSys)
207 
208  # Build the Shapelet objects
209  ctrl = config.makeControl()
210  baseName = "modelfit_CModel"
211  nonlinearKeys = ["{}_{{model}}_nonlinear_{p}".format(baseName, p=p)
212  for p in range(3)]
213  fixedKeys = ["{}_{{model}}_fixed_{p}".format(baseName, p=p)
214  for p in range(2)]
215  fluxKey = "{}_{{model}}_instFlux".format(baseName)
216 
217  # fetch the aperture corrections, if this fails set it to one
218  try:
219  apCorr = record.get("{}_apCorr".format(baseName))
220  except Exception:
221  print("Warning, problem retrieving aperture correction, using a value"
222  " of 1")
223  apCorr = 1
224 
225  # Get parameters for the dev model
226  devNonLinearParams = np.array([record.get(key.format(model="dev"))
227  for key in nonlinearKeys])
228  devFixedParams = np.array([record.get(key.format(model="dev"))
229  for key in fixedKeys])
230  devAmp = np.array([record.get(fluxKey.format(model="dev"))])
231  devAmp /= apCorr
232  devShapelet = ctrl.dev.getModel().makeShapeletFunction(devNonLinearParams,
233  devAmp,
234  devFixedParams)
235  devShapelet.transformInPlace(fitSysToMeasSys.geometric)
236 
237  # Get parameters for the exp model
238  expNonLinearParams = np.array([record.get(key.format(model="exp"))
239  for key in nonlinearKeys])
240  expFixedParams = np.array([record.get(key.format(model="exp"))
241  for key in fixedKeys])
242  expAmp = np.array([record.get(fluxKey.format(model="exp"))])
243  expAmp /= apCorr
244  expShapelet = ctrl.exp.getModel().makeShapeletFunction(expNonLinearParams,
245  expAmp,
246  expFixedParams)
247  expShapelet.transformInPlace(fitSysToMeasSys.geometric)
248 
249  # Get joint shapelet model
250  fracDev = record.get("{}_fracDev".format(baseName))
251  jointFlux = np.array([record.get("{}_instFlux".format(baseName))])
252  jointFlux /= apCorr
253  devJointShapelet = ctrl.dev.getModel()\
254  .makeShapeletFunction(devNonLinearParams, jointFlux*fracDev,
255  devFixedParams)
256  devJointShapelet.transformInPlace(fitSysToMeasSys.geometric)
257 
258  expJointShapelet = ctrl.exp.getModel()\
259  .makeShapeletFunction(expNonLinearParams, jointFlux*(1-fracDev),
260  expFixedParams)
261  expJointShapelet.transformInPlace(fitSysToMeasSys.geometric)
262 
263  return devShapelet, expShapelet, devJointShapelet, expJointShapelet
int min
int max
Class that maps MultiShapeletFunction objects to fields in afw::table objects.
Definition: FunctorKeys.h:161
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
Point< double, 2 > Point2D
Definition: Point.h:324
void addToImage(std::shared_ptr< afw::image::Image< double >> image, std::vector< std::shared_ptr< afw::image::Image< double >>> const &imgVector, std::vector< double > const &weightVector)
Definition: CoaddPsf.cc:201
def buildCModelImages(exposure, record, config)
def reconstructCModel(exposure, record, config)
def displayReconstructedCmodel(exposure, record, config, display="mpl")
def displayReconstrucedCmodelMpl(exposure, record, config)
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174