Coverage for python / lsst / analysis / tools / actions / plot / multiVisitCoveragePlot.py: 8%

1# This file is part of analysis_tools. 

2# 

3# Developed for the LSST Data Management System. 

4# This product includes software developed by the LSST Project 

5# (https://www.lsst.org). 

6# See the COPYRIGHT file at the top-level directory of this distribution 

7# for details of code ownership. 

8# 

9# This program is free software: you can redistribute it and/or modify 

10# it under the terms of the GNU General Public License as published by 

11# the Free Software Foundation, either version 3 of the License, or 

12# (at your option) any later version. 

13# 

14# This program is distributed in the hope that it will be useful, 

15# but WITHOUT ANY WARRANTY; without even the implied warranty of 

16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

17# GNU General Public License for more details. 

18# 

19# You should have received a copy of the GNU General Public License 

20# along with this program. If not, see <https://www.gnu.org/licenses/>. 

21from __future__ import annotations 

22 

23__all__ = ("MultiVisitCoveragePlot",) 

24 

25import logging 

26from collections.abc import Mapping 

27from typing import cast 

28 

29import matplotlib as mpl 

30import matplotlib.pyplot as plt 

31import numpy as np 

32import pandas as pd 

33from matplotlib.figure import Figure 

34from matplotlib.ticker import FormatStrFormatter 

35 

36from lsst.afw.cameraGeom import FOCAL_PLANE, Camera, DetectorType 

37from lsst.pex.config import Config, DictField, Field, ListField 

38from lsst.skymap import BaseSkyMap 

39 

40from ...interfaces import KeyedData, KeyedDataSchema, PlotAction, Scalar, Vector 

41from ...math import nanMax, nanMean, nanMedian, nanMin 

42from ..vector.selectors import RangeSelector 

43from .plotUtils import mkColormap, plotProjectionWithBinning 

44 

45log = logging.getLogger(__name__) 

46 

47 

48class MultiVisitCoveragePlot(PlotAction): 

49 """Plot the coverage for a set of visits.""" 

50 

51 plotName = Field[str]( 

52 doc="The name for the plot.", 

53 optional=False, 

54 ) 

55 cmapName = Field[str]( 

56 doc="Name of the color map to be used if not using the default color-blind friendly " 

57 "orange/blue default (used if this is set to `None`). Any name available via " 

58 "`matplotlib.cm` may be used.", 

59 default=None, 

60 optional=True, 

61 ) 

62 projection = Field[str]( 

63 doc='Projection to plot. Currently only "raDec" and "focalPlane" are permitted. ' 

64 "In either case, one point is plotted per visit/detector combination.", 

65 default="raDec", 

66 ) 

67 nBins = Field[int]( 

68 doc="Number of bins to use within the effective plot ranges along the spatial directions. " 

69 'Only used in the "raDec" projection (for the "focalPlane" projection, the binning is ' 

70 "effectively one per detector).", 

71 default=25, 

72 ) 

73 nPointBinThresh = Field[int]( 

74 doc="Threshold number of data points above which binning of the data will be performed in " 

75 'the RA/Dec projection. If ``projection`` is "focalPlane", the per-detector nPoint ' 

76 "threshold is nPointMinThresh/number of science detectors in the given ``camera``.", 

77 default=400, 

78 ) 

79 showExtremeOutliers = Field[bool]( 

80 doc="Use overlaid scatter points to show the x-y positions of the 15% most extreme values?", 

81 default=True, 

82 ) 

83 unitsDict = DictField[str, str]( 

84 doc="A dict mapping a parameter to its appropriate units (for label plotting).", 

85 default={ 

86 "astromOffsetMean": "arcsec", 

87 "astromOffsetStd": "arcsec", 

88 "psfSigma": "pixels", 

89 "skyBg": "counts", 

90 "skyNoise": "counts", 

91 "visit": "number", 

92 "detector": "number", 

93 "zenithDistance": "deg", 

94 "zeroPoint": "mag", 

95 "ra": "deg", 

96 "decl": "deg", 

97 "xFp": "mm", 

98 "yFp": "mm", 

99 "medianE": "", 

100 "psfStarScaledDeltaSizeScatter": "", 

101 "psfTraceRadiusDelta": "pixel", 

102 }, 

103 ) 

104 sortedFullBandList = ListField[str]( 

105 doc="List of all bands that could, in principle, but do not have to, exist in data table. " 

106 "The sorting of the plot panels will follow this list (typically by wavelength).", 

107 default=["u", "g", "r", "i", "z", "y", "N921"], 

108 ) 

109 bandLabelColorDict = DictField[str, str]( 

110 doc="A dict mapping which color to use for the labels of a given band.", 

111 default={ 

112 "u": "tab:purple", 

113 "g": "tab:blue", 

114 "r": "tab:green", 

115 "i": "gold", 

116 "z": "tab:orange", 

117 "y": "tab:red", 

118 "N921": "tab:pink", 

119 }, 

120 ) 

121 vectorsToLoadList = ListField[str]( 

122 doc="List of columns to load from input table.", 

123 default=[ 

124 "visitId", 

125 "detector", 

126 "band", 

127 "ra", 

128 "decl", 

129 "zeroPoint", 

130 "psfSigma", 

131 "skyBg", 

132 "astromOffsetMean", 

133 "psfStarDeltaE1Median", 

134 "psfStarDeltaE2Median", 

135 "psfStarScaledDeltaSizeScatter", 

136 "psfTraceRadiusDelta", 

137 "llcra", 

138 "lrcra", 

139 "ulcra", 

140 "urcra", 

141 "llcdec", 

142 "lrcdec", 

143 "ulcdec", 

144 "urcdec", 

145 "xSize", 

146 "ySize", 

147 ], 

148 ) 

149 parametersToPlotList = ListField[str]( 

150 doc="List of parameters to plot. They are plotted along rows and the columns " 

151 "plot these parameters for each band.", 

152 default=[ 

153 "psfSigma", 

154 "astromOffsetMean", 

155 "medianE", 

156 "psfStarScaledDeltaSizeScatter", 

157 "skyBg", 

158 "zeroPoint", 

159 ], 

160 ) 

161 tractsToPlotList = ListField[int]( 

162 doc="List of tracts within which to limit the RA and Dec limits of the plot.", 

163 default=None, 

164 optional=True, 

165 ) 

166 trimToTract = Field[bool]( 

167 doc="Whether to trim the plot limits to the tract limit(s). Otherwise, plot " 

168 "will be trimmed to data limits (both will be expanded in the smaller range " 

169 "direction for an equal aspect square plot).", 

170 default=False, 

171 ) 

172 doScatterInRaDec = Field[bool]( 

173 doc="Whether to scatter the points in RA/Dec before plotting. This may be useful " 

174 "for visualization for surveys with tight dithering patterns.", 

175 default=False, 

176 ) 

177 plotAllTractOutlines = Field[bool]( 

178 doc="Whether to plot tract outlines for all tracts within the plot limits " 

179 "(regardless if they have any data in them).", 

180 default=False, 

181 ) 

182 raDecLimitsDict = DictField[str, float]( 

183 doc="A dict mapping the RA/Dec limits to apply to the plot. Set to `None` to use " 

184 "base limits on the default or the other config options. The dict must contain " 

185 "the keys raMin, raMax, decMin, decMax, e.g. " 

186 'raDecLimitsDict = {"raMin": 0, "raMax": 360, "decMin": -90, "decMax": 90}. ' 

187 "Not compatible with ``trimToTract`` or ``tractsToPlotList`` (i.e. the latter two " 

188 "will be ignored if the dict is not `None`).", 

189 default=None, 

190 optional=True, 

191 ) 

192 plotDetectorOutline = Field[bool]( 

193 doc="Whether to plot a shaded outline of the detector size in the RA/Dec projection" 

194 "for reference. Ignored if ``projection`` is not raDec or no camera is provided " 

195 "in the inputs.", 

196 default=False, 

197 ) 

198 

199 def getInputSchema(self) -> KeyedDataSchema: 

200 base: list[tuple[str, type[Vector] | type[Scalar]]] = [] 

201 for vector in self.vectorsToLoadList: 

202 base.append((vector, Vector)) 

203 return base 

204 

205 def __call__(self, data: KeyedData, **kwargs) -> Mapping[str, Figure] | Figure: 

206 self._validateInput(data, **kwargs) 

207 return self.makePlot(data, **kwargs) 

208 

209 def _validateInput(self, data: KeyedData, **kwargs) -> None: 

210 """NOTE currently can only check that something is not a Scalar, not 

211 check that the data is consistent with Vector. 

212 """ 

213 needed = self.getFormattedInputSchema(**kwargs) 

214 if remainder := {key.format(**kwargs) for key, _ in needed} - { 

215 key.format(**kwargs) for key in data.keys() 

216 }: 

217 raise ValueError(f"Task needs keys {remainder} but they were not found in input.") 

218 for name, typ in needed: 

219 isScalar = issubclass((colType := type(data[name.format(**kwargs)])), Scalar) 

220 if isScalar and typ != Scalar: 

221 raise ValueError(f"Data keyed by {name} has type {colType} but action requires type {typ}.") 

222 

223 if "medianE" in self.parametersToPlotList: 

224 if not all( 

225 vector in self.vectorsToLoadList 

226 for vector in ["psfStarDeltaE1Median", "psfStarDeltaE2Median"] 

227 ): 

228 raise RuntimeError( 

229 "If medianE is to be plotted, both psfStarDeltaE1Median and " 

230 "psfStarDeltaE2Median must be included in vectorsToLoadList." 

231 ) 

232 if self.raDecLimitsDict is not None: 

233 requiredKeys = set(["raMin", "raMax", "decMin", "decMax"]) 

234 if requiredKeys != set(self.raDecLimitsDict.keys()): 

235 raise RuntimeError( 

236 f"The following keys (and only these) are required in raDecLimitDict: {requiredKeys}." 

237 f"The dict provided gave {set(self.raDecLimitsDict.keys())}." 

238 ) 

239 

240 def makePlot( 

241 self, 

242 data: KeyedData, 

243 plotInfo: Mapping[str, str] | None = None, 

244 camera: Camera | None = None, 

245 skymap: BaseSkyMap | None = None, 

246 calibrateConfig: Config | None = None, 

247 makeWarpConfig: Config | None = None, 

248 **kwargs, 

249 ) -> Figure: 

250 """Make an Nband x Nparameter panel multi-visit coverage plot. 

251 

252 The panels rows are for different bands,sorted according to the order 

253 in config ``sortedFullBandList``. The columns are per-parameter, 

254 plotted in the order given by the config ``parametersToPlotList``. 

255 Red "over" colors indicate thresholds in play in the data processing 

256 (e.g. used to select against data of sufficient quality). 

257 

258 Parameters 

259 ---------- 

260 data : `lsst.analysis.tools.interfaces.KeyedData` 

261 The key-based catalog of data to plot. 

262 plotInfo : `dict` [`str`], optional 

263 A dictionary of information about the data being plotted with (at 

264 least) keys: 

265 

266 `"run"` 

267 Output run for the plots (`str`). 

268 `"tableName"` 

269 Name of the table from which results are taken (`str`). 

270 camera : `lsst.afw.cameraGeom.Camera`, optional 

271 The camera object associated with the data. This is to enable the 

272 conversion of to focal plane coordinates (if needed, i.e. for the 

273 focal plane projection version of this plot) and to obtain the 

274 projected (in RA/Dec) size of a typical detector (for reference in 

275 the raDec projection plots when requested, i.e. if the config 

276 ``plotDetectorOutline`` is `True`). 

277 skymap : `lsst.skymap.BaseSkyMap`, optional 

278 The sky map used for this dataset. If a specific tract[List] is 

279 provided, this is used to determine the appropriate RA & Dec limits 

280 to downselect the data to within those limits for plotting. 

281 calibrateConfig : `lsst.pex.config.Config`, optional 

282 The persisted config used in the calibration task for the given 

283 collection. Used to introspect threshold values used in the run. 

284 makeWarpConfig : `lsst.pex.config.Config`, optional 

285 The persisted config used in the makeDirectWarp task 

286 for the given collection. Used to introspect threshold values used 

287 in the run. 

288 

289 Returns 

290 ------- 

291 fig : `matplotlib.figure.Figure` 

292 The resulting figure. 

293 """ 

294 mpl.rcParams["figure.dpi"] = 350 

295 mpl.rcParams["font.size"] = 5 

296 mpl.rcParams["xtick.labelsize"] = 5 

297 mpl.rcParams["ytick.labelsize"] = 5 

298 mpl.rcParams["xtick.major.size"] = 1.5 

299 mpl.rcParams["ytick.major.size"] = 1.5 

300 mpl.rcParams["xtick.major.width"] = 0.5 

301 mpl.rcParams["ytick.major.width"] = 0.5 

302 

303 if "medianE" in self.parametersToPlotList: 

304 data["medianE"] = np.sqrt( 

305 data["psfStarDeltaE1Median"] ** 2.0 + data["psfStarDeltaE2Median"] ** 2.0 

306 ) 

307 

308 # Make sure all columns requested actually exist in ccdVisitTable. 

309 notFoundKeys = [] 

310 for zKey in self.parametersToPlotList: 

311 if zKey not in data.keys(): 

312 notFoundKeys.append(zKey) 

313 if len(notFoundKeys) > 0: 

314 raise RuntimeError(f"Requested column(s) {notFoundKeys} is(are) not in the data table.") 

315 

316 maxMeanDistanceArcsec = calibrateConfig.astrometry.maxMeanDistanceArcsec # type: ignore 

317 

318 if makeWarpConfig is None: 

319 maxEllipResidual = 0.007 

320 maxScaledSizeScatter = 0.019 

321 maxPsfTraceRadiusDelta = 0.7 

322 else: 

323 maxEllipResidual = makeWarpConfig.select.value.maxEllipResidual # type: ignore 

324 maxScaledSizeScatter = makeWarpConfig.select.value.maxScaledSizeScatter # type: ignore 

325 maxPsfTraceRadiusDelta = makeWarpConfig.select.value.maxPsfTraceRadiusDelta # type: ignore 

326 

327 cameraName = "" if camera is None else camera.getName() 

328 if self.projection == "focalPlane": 

329 if camera is None: 

330 raise RuntimeError("Must have an input camera if plotting focalPlane projection.") 

331 xKey = "xFp" 

332 yKey = "yFp" 

333 # Add the detector center in Focal Plane coords to the data table. 

334 detFpDict = {} 

335 for det in camera: # type: ignore 

336 if det.getType() == DetectorType.SCIENCE: 

337 detFpDict[det.getId()] = det.getCenter(FOCAL_PLANE) 

338 log.info(f"Number of SCIENCE detectors in {cameraName} camera = {len(detFpDict)}") 

339 xFpList = [] 

340 yFpList = [] 

341 for det in data["detector"]: # type: ignore 

342 xFpList.append(detFpDict[det].getX()) 

343 yFpList.append(detFpDict[det].getY()) 

344 

345 data["xFp"] = xFpList # type: ignore 

346 data["yFp"] = yFpList # type: ignore 

347 

348 corners = camera[0].getCorners(FOCAL_PLANE) # type: ignore 

349 xCorners, yCorners = zip(*corners) 

350 xScatLen = 0.4 * (nanMax(xCorners) - nanMin(xCorners)) 

351 yScatLen = 0.4 * (nanMax(yCorners) - nanMin(yCorners)) 

352 tractList: list[int] = [] 

353 elif self.projection == "raDec": 

354 xKey = "ra" 

355 yKey = "decl" 

356 xScatLen, yScatLen = 0, 0 

357 # Use downselector without limits to get rid of any non-finite 

358 # RA/Dec entries. 

359 data = self._planeAreaSelector(data, xKey=xKey, yKey=yKey) 

360 if self.tractsToPlotList is not None and len(self.tractsToPlotList) > 0: 

361 tractList = self.tractsToPlotList 

362 else: 

363 ras = data["ra"] 

364 decs = data["decl"] 

365 tractList = list(set(skymap.findTractIdArray(ras, decs, degrees=True))) # type: ignore 

366 log.info(f"List of tracts overlapping data: {tractList}") 

367 tractLimitsDict = self._getTractLimitsDict(skymap, tractList) 

368 else: 

369 raise ValueError(f"Unknown projection: {self.projection}") 

370 

371 perTract = True if self.tractsToPlotList is not None and self.projection == "raDec" else False 

372 

373 if perTract and len(self.tractsToPlotList) > 0: 

374 data = self._trimDataToTracts(data, skymap) 

375 nData = len(cast(Vector, data[list(data.keys())[0]])) 

376 if nData == 0: 

377 raise RuntimeError( 

378 f"No data to plot. Did your tract selection of {self.tractsToPlotList} remove all data?" 

379 ) 

380 

381 if self.doScatterInRaDec: 

382 raRange = max(cast(Vector, data["ra"])) - min(cast(Vector, data["ra"])) 

383 decRange = max(cast(Vector, data["decl"])) - min(cast(Vector, data["decl"])) 

384 scatRad = max(0.05 * max(raRange, decRange), 0.12) # min is of order of an LSSTCam detector 

385 log.info(f"Scattering data in RA/Dec within radius {scatRad:.3f} (deg)") 

386 

387 dataDf = pd.DataFrame(data) 

388 nDataId = len(dataDf) 

389 log.info("Number of detector data points: %i", nDataId) 

390 if nDataId == 0: 

391 raise RuntimeError("No data to plot. Exiting...") 

392 nVisit = len(set(dataDf["visitId"])) 

393 

394 # Make a sorted list of the bands that exist in the data table. 

395 dataBandList = list(set(dataDf["band"])) 

396 bandList = [band for band in self.sortedFullBandList if band in dataBandList] 

397 missingBandList = list(set(dataBandList) - set(self.sortedFullBandList)) 

398 if len(missingBandList) > 0: 

399 log.warning( 

400 f"The band(s) {missingBandList} are not included in self.sortedFullBandList. Please add them " 

401 "so they get sorted properly (namely by wavelength). For now, they will just be " 

402 "appended to the end of the list of those that could be sorted. You may also " 

403 "wish to give them entries in self.bandLabelColorList to specify the label color " 

404 "(otherwise, if not present, it defaults to teal)." 

405 ) 

406 bandList.extend(missingBandList) 

407 log.info(f"Sorted list of existing bands: {bandList}") 

408 

409 ptSize = min(5, max(0.3, 600.0 / ((nDataId / len(bandList)) ** 0.5))) 

410 if self.doScatterInRaDec: 

411 ptSize = min(3, 0.7 * ptSize) 

412 

413 nRow = len(bandList) 

414 nCol = len(self.parametersToPlotList) 

415 colMultiplier = 4 if nCol == 1 else 2.5 

416 fig, axes = plt.subplots( 

417 nRow, 

418 nCol, 

419 figsize=(int(colMultiplier * nCol), 2 * nRow), 

420 constrained_layout=True, 

421 ) 

422 

423 # Select reasonable limits for colormaps so they can be matched for 

424 # all bands. 

425 nPercent = max(2, int(0.02 * nDataId)) 

426 vMinDict = {} 

427 vMaxDict = {} 

428 for zKey in self.parametersToPlotList: 

429 zKeySorted = dataDf[zKey].sort_values() 

430 zKeySorted = zKeySorted[np.isfinite(zKeySorted)] 

431 vMinDict[zKey] = nanMean(zKeySorted.head(nPercent)) 

432 if zKey == "medianE": 

433 vMaxDict[zKey] = maxEllipResidual 

434 elif zKey == "psfStarScaledDeltaSizeScatter": 

435 vMaxDict[zKey] = maxScaledSizeScatter 

436 elif zKey == "astromOffsetMean" and self.projection != "raDec": 

437 vMaxDict[zKey] = min(maxMeanDistanceArcsec, 1.1 * nanMean(zKeySorted.tail(nPercent))) 

438 elif zKey == "maxPsfTraceRadiusDelta": 

439 vMaxDict[zKey] = maxPsfTraceRadiusDelta 

440 else: 

441 vMaxDict[zKey] = nanMean(zKeySorted.tail(nPercent)) 

442 

443 for iRow, band in enumerate(bandList): 

444 dataBand = dataDf[dataDf["band"] == band].copy() 

445 

446 nDataIdBand = len(dataBand) 

447 nVisitBand = len(set(dataBand["visitId"])) 

448 if nDataIdBand < 2: 

449 log.warning(f"Fewer than 2 points to plot for {band}. Skipping...") 

450 continue 

451 

452 for zKey in self.parametersToPlotList: 

453 # Don't match the colorbars when it doesn't make sense to. 

454 if zKey in ["skyBg", "zeroPoint"]: 

455 nPercent = max(2, int(0.02 * nDataIdBand)) 

456 zKeySorted = dataBand[zKey].sort_values() 

457 zKeySorted = zKeySorted[np.isfinite(zKeySorted)] 

458 vMinDict[zKey] = nanMean(zKeySorted.head(nPercent)) 

459 vMaxDict[zKey] = nanMean(zKeySorted.tail(nPercent)) 

460 

461 # Scatter the plots within the detector for focal plane plots. 

462 if self.doScatterInRaDec: 

463 scatRads = scatRad * np.sqrt(np.random.uniform(size=nDataIdBand)) 

464 scatTheta = 2.0 * np.pi * np.random.uniform(size=nDataIdBand) 

465 xScatter = scatRads * np.cos(scatTheta) 

466 yScatter = scatRads * np.sin(scatTheta) 

467 xLabel = xKey + " + rand(scatter)" 

468 yLabel = yKey + " + rand(scatter)" 

469 else: 

470 xScatter = np.random.uniform(-xScatLen, xScatLen, len(dataBand[xKey])) 

471 yScatter = np.random.uniform(-yScatLen, yScatLen, len(dataBand[yKey])) 

472 xLabel = xKey 

473 yLabel = yKey 

474 dataBand["xScat"] = dataBand[xKey] + xScatter 

475 dataBand["yScat"] = dataBand[yKey] + yScatter 

476 # Accommodate odd number of quarter-turn rotated detectors. 

477 if self.projection == "focalPlane": 

478 for det in camera: # type: ignore 

479 if det.getOrientation().getNQuarter() % 2 != 0: 

480 detId = int(det.getId()) 

481 xFpRot = dataBand.loc[dataBand.detector == detId, xKey] 

482 yFpRot = dataBand.loc[dataBand.detector == detId, yKey] 

483 xScatRot = dataBand.loc[dataBand.detector == detId, "xScat"] 

484 yScatRot = dataBand.loc[dataBand.detector == detId, "yScat"] 

485 dataBand.loc[dataBand.detector == detId, "xScat"] = xFpRot + (yScatRot - yFpRot) 

486 dataBand.loc[dataBand.detector == detId, "yScat"] = yFpRot + (xScatRot - xFpRot) 

487 

488 if band not in self.bandLabelColorDict: 

489 log.warning( 

490 f"The band {band} is not included in the bandLabelColorList config. Please add it " 

491 "to specify the label color (otherwise, it defaults to teal)." 

492 ) 

493 color = self.bandLabelColorDict[band] if band in self.bandLabelColorDict else "teal" 

494 fontDict = {"fontsize": 5, "color": color} 

495 

496 for iCol, zKey in enumerate(self.parametersToPlotList): 

497 if self.cmapName is None: 

498 cmap = mkColormap(["darkOrange", "thistle", "midnightblue"]) 

499 else: 

500 cmap = mpl.cm.get_cmap(self.cmapName).copy() 

501 

502 if zKey in ["medianE", "psfStarScaledDeltaSizeScatter", "psfTraceRadiusDelta"]: 

503 cmap.set_over("red") 

504 elif ( 

505 zKey in ["astromOffsetMean"] 

506 and self.projection != "raDec" 

507 and vMaxDict[zKey] >= maxMeanDistanceArcsec 

508 ): 

509 cmap.set_over("red") 

510 else: 

511 if self.cmapName is None: 

512 cmap.set_over("black") 

513 else: 

514 cmap.set_over("lemonchiffon") 

515 

516 if zKey in ["psfSigma"]: 

517 cmap.set_under("red") 

518 else: 

519 if self.cmapName is None: 

520 cmap.set_under("lemonchiffon") 

521 else: 

522 cmap.set_under("black") 

523 

524 titleStr = f"band: {band} nVisit: {nVisitBand} nData: {nDataIdBand}" 

525 

526 ax = axes[iRow, iCol] if axes.ndim > 1 else axes[max(iRow, iCol)] 

527 ax.set_title(f"{titleStr}", loc="left", fontdict=fontDict, pad=2) 

528 ax.set_xlabel(f"{xLabel} ({self.unitsDict[xKey]})", labelpad=0) 

529 ax.set_ylabel(f"{yLabel} ({self.unitsDict[yKey]})", labelpad=1) 

530 ax.set_aspect("equal") 

531 ax.tick_params("x", labelrotation=45, pad=0) 

532 

533 if self.projection == "focalPlane": 

534 for det in camera: # type: ignore 

535 if det.getType() == DetectorType.SCIENCE: 

536 corners = det.getCorners(FOCAL_PLANE) 

537 xCorners, yCorners = zip(*corners) 

538 xFpMin, xFpMax = min(xCorners), max(xCorners) 

539 yFpMin, yFpMax = min(yCorners), max(yCorners) 

540 detId = int(det.getId()) 

541 perDetData = dataBand[dataBand["detector"] == detId] 

542 if len(perDetData) < 1: 

543 log.debug(f"No data to plot for detector {detId}. Skipping...") 

544 continue 

545 if sum(np.isfinite(perDetData[zKey])) < 1: 

546 log.debug(f"No finited data to plot for detector {detId}. Skipping...") 

547 continue 

548 pcm = plotProjectionWithBinning( 

549 ax, 

550 perDetData["xScat"].values, 

551 perDetData["yScat"].values, 

552 perDetData[zKey].values, 

553 cmap, 

554 xFpMin, 

555 xFpMax, 

556 yFpMin, 

557 yFpMax, 

558 xNumBins=1, 

559 fixAroundZero=False, 

560 nPointBinThresh=max(1, int(self.nPointBinThresh / len(detFpDict))), 

561 isSorted=False, 

562 vmin=vMinDict[zKey], 

563 vmax=vMaxDict[zKey], 

564 showExtremeOutliers=self.showExtremeOutliers, 

565 scatPtSize=ptSize, 

566 ) 

567 if self.projection == "raDec": 

568 raMin, raMax = min(cast(Vector, data["ra"])), max(cast(Vector, data["ra"])) 

569 decMin, decMax = min(cast(Vector, data["decl"])), max(cast(Vector, data["decl"])) 

570 raMin, raMax, decMin, decMax = _setLimitsToEqualRatio(raMin, raMax, decMin, decMax) 

571 pcm = plotProjectionWithBinning( 

572 ax, 

573 dataBand["xScat"].values, 

574 dataBand["yScat"].values, 

575 dataBand[zKey].values, 

576 cmap, 

577 raMin, 

578 raMax, 

579 decMin, 

580 decMax, 

581 xNumBins=self.nBins, 

582 fixAroundZero=False, 

583 nPointBinThresh=self.nPointBinThresh, 

584 isSorted=False, 

585 vmin=vMinDict[zKey], 

586 vmax=vMaxDict[zKey], 

587 showExtremeOutliers=self.showExtremeOutliers, 

588 scatPtSize=ptSize, 

589 ) 

590 

591 # Make sure all panels get the same axis limits and always make 

592 # equal axis ratio panels. 

593 if iRow == 0 and iCol == 0: 

594 if self.raDecLimitsDict is not None and self.projection == "raDec": 

595 xLimMin = self.raDecLimitsDict["raMin"] # type: ignore 

596 xLimMax = self.raDecLimitsDict["raMax"] # type: ignore 

597 yLimMin = self.raDecLimitsDict["decMin"] # type: ignore 

598 yLimMax = self.raDecLimitsDict["decMax"] # type: ignore 

599 else: 

600 xLimMin, xLimMax = ax.get_xlim() 

601 yLimMin, yLimMax = ax.get_ylim() 

602 if self.projection == "focalPlane": 

603 minDim = ( 

604 max( 

605 camera.getFpBBox().getWidth(), # type: ignore 

606 camera.getFpBBox().getHeight(), # type: ignore 

607 ) 

608 / 2 

609 ) 

610 xLimMin = min(-minDim, xLimMin) 

611 xLimMax = max(minDim, xLimMax) 

612 if self.trimToTract: 

613 for tract, tractLimits in tractLimitsDict.items(): 

614 xLimMin = min(xLimMin, min(tractLimits["ras"])) 

615 xLimMax = max(xLimMax, max(tractLimits["ras"])) 

616 yLimMin = min(yLimMin, min(tractLimits["decs"])) 

617 yLimMax = max(yLimMax, max(tractLimits["decs"])) 

618 xDelta = xLimMax - xLimMin 

619 xLimMin -= 0.04 * xDelta 

620 xLimMax += 0.04 * xDelta 

621 

622 xLimMin, xLimMax, yLimMin, yLimMax = _setLimitsToEqualRatio( 

623 xLimMin, xLimMax, yLimMin, yLimMax 

624 ) 

625 limRange = xLimMax - xLimMin 

626 yTickFmt = _tickFormatter(yLimMin * 10) 

627 

628 if self.plotAllTractOutlines and self.projection == "raDec": 

629 allTractsList = [] 

630 for tractInfo in skymap: # type: ignore 

631 centerRa = tractInfo.getCtrCoord()[0].asDegrees() 

632 centerDec = tractInfo.getCtrCoord()[1].asDegrees() 

633 if ( 

634 centerRa > xLimMin 

635 and centerRa < xLimMax 

636 and centerDec > yLimMin 

637 and centerDec < yLimMax 

638 ): 

639 allTractsList.append(int(tractInfo.getId())) 

640 tractLimitsDict = self._getTractLimitsDict(skymap, allTractsList) 

641 

642 upperHandles = [] 

643 if self.doScatterInRaDec and self.projection == "raDec": 

644 patch = mpl.patches.Circle( 

645 (xLimMax - 1.5 * scatRad, yLimMax - 1.5 * scatRad), 

646 radius=scatRad, 

647 facecolor="gray", 

648 edgecolor="None", 

649 alpha=0.2, 

650 label="scatter area", 

651 ) 

652 ax.add_patch(patch) 

653 upperHandles.append(patch) 

654 

655 # Add a shaded area of the size of a detector for reference. 

656 if self.plotDetectorOutline and self.projection == "raDec": 

657 if camera is None: 

658 if iRow == 0 and iCol == 0: 

659 log.warning( 

660 "Config plotDetectorOutline is True, but no camera was provided. " 

661 "Reference detector outline will not be included in the plot." 

662 ) 

663 else: 

664 # Calculate area of polygon with known vertices. 

665 x1, x2, x3, x4 = ( 

666 dataBand["llcra"], 

667 dataBand["lrcra"], 

668 dataBand["urcra"], 

669 dataBand["ulcra"], 

670 ) 

671 y1, y2, y3, y4 = ( 

672 dataBand["llcdec"], 

673 dataBand["lrcdec"], 

674 dataBand["urcdec"], 

675 dataBand["ulcdec"], 

676 ) 

677 areaDeg = ( 

678 np.abs( 

679 (x1 * y2 - y1 * x2) 

680 + (x2 * y3 - y2 * x3) 

681 + (x3 * y4 - y3 * x4) 

682 + (x4 * y1 - y4 * x1) 

683 ) 

684 / 2.0 

685 ) 

686 detScaleDeg = np.sqrt(areaDeg / (dataBand["xSize"] * dataBand["ySize"])) 

687 detWidthDeg = nanMedian(detScaleDeg * dataBand["xSize"]) 

688 detHeightDeg = nanMedian(detScaleDeg * dataBand["ySize"]) 

689 

690 patch = mpl.patches.Rectangle( 

691 (xLimMax - 0.02 * limRange - detWidthDeg, yLimMin + 0.03 * limRange), 

692 detWidthDeg, 

693 detHeightDeg, 

694 facecolor="turquoise", 

695 alpha=0.3, 

696 label="det size", 

697 ) 

698 ax.add_patch(patch) 

699 upperHandles.append(patch) 

700 

701 if self.projection == "raDec": 

702 if iRow == 0 and iCol == 0 and len(upperHandles) > 0: 

703 ax.legend( 

704 handles=upperHandles, 

705 loc="upper left", 

706 bbox_to_anchor=(0, 1.17 + 0.05 * len(upperHandles)), 

707 edgecolor="black", 

708 framealpha=0.4, 

709 fontsize=5, 

710 ) 

711 

712 # Overplot tract outlines if tracts were specified, but only if 

713 # the plot limits span fewer than the width of 5 tracts 

714 # (otherwise the labels will be too crowded to be useful). 

715 if len(tractList) > 0: 

716 tractInfo = skymap[tractList[0]] # type: ignore 

717 tractWidthDeg = tractInfo.outer_sky_polygon.getBoundingBox().getWidth().asDegrees() 

718 if limRange <= 5 * tractWidthDeg: 

719 deltaLim = 0.05 * limRange 

720 for tract, tractLimits in tractLimitsDict.items(): 

721 centerRa = tractLimits["center"][0] 

722 centerDec = tractLimits["center"][1] 

723 if ( 

724 centerRa > xLimMin + deltaLim 

725 and centerRa < xLimMax - deltaLim 

726 and centerDec > yLimMin + deltaLim 

727 and centerDec < yLimMax - deltaLim 

728 ): 

729 ax.plot(tractLimits["ras"], tractLimits["decs"], color="dimgray", lw=0.4) 

730 fontSize = 3 if limRange < 20 else 2 

731 ax.annotate( 

732 str(tract), 

733 tractLimits["center"], 

734 va="center", 

735 ha="center", 

736 color="dimgray", 

737 fontsize=fontSize, 

738 annotation_clip=True, 

739 path_effects=[mpl.patheffects.withStroke(linewidth=1, foreground="w")], 

740 ) 

741 if self.projection == "raDec": 

742 ax.set_xlim(xLimMax, xLimMin) 

743 else: 

744 ax.set_xlim(xLimMin, xLimMax) 

745 ax.set_ylim(yLimMin, yLimMax) 

746 ax.yaxis.set_major_formatter(yTickFmt) 

747 

748 # Get a tick formatter that will give all anticipated value 

749 # ranges the same length. This is so that the label padding 

750 # has the same effect on all colorbars. 

751 value = vMaxDict[zKey] if np.abs(vMaxDict[zKey]) > np.abs(vMinDict[zKey]) else vMinDict[zKey]