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 typing import List, Mapping, Optional, cast 

27 

28import matplotlib as mpl 

29import matplotlib.pyplot as plt 

30import numpy as np 

31import pandas as pd 

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

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

34from lsst.skymap import BaseSkyMap 

35from matplotlib.figure import Figure 

36from matplotlib.ticker import FormatStrFormatter 

37 

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

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

40from ..vector.selectors import RangeSelector 

41from .plotUtils import mkColormap, plotProjectionWithBinning 

42 

43log = logging.getLogger(__name__) 

44 

45 

46class MultiVisitCoveragePlot(PlotAction): 

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

48 

49 plotName = Field[str]( 

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

51 optional=False, 

52 ) 

53 cmapName = Field[str]( 

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

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

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

57 default=None, 

58 optional=True, 

59 ) 

60 projection = Field[str]( 

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

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

63 default="raDec", 

64 ) 

65 nBins = Field[int]( 

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

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

68 "effectively one per detector).", 

69 default=25, 

70 ) 

71 nPointBinThresh = Field[int]( 

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

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

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

75 default=400, 

76 ) 

77 showExtremeOutliers = Field[bool]( 

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

79 default=True, 

80 ) 

81 unitsDict = DictField[str, str]( 

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

83 default={ 

84 "astromOffsetMean": "arcsec", 

85 "astromOffsetStd": "arcsec", 

86 "psfSigma": "pixels", 

87 "skyBg": "counts", 

88 "skyNoise": "counts", 

89 "visit": "number", 

90 "detector": "number", 

91 "zenithDistance": "deg", 

92 "zeroPoint": "mag", 

93 "ra": "deg", 

94 "decl": "deg", 

95 "xFp": "mm", 

96 "yFp": "mm", 

97 "medianE": "", 

98 "psfStarScaledDeltaSizeScatter": "", 

99 "psfTraceRadiusDelta": "pixel", 

100 }, 

101 ) 

102 sortedFullBandList = ListField[str]( 

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

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

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

106 ) 

107 bandLabelColorDict = DictField[str, str]( 

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

109 default={ 

110 "u": "tab:purple", 

111 "g": "tab:blue", 

112 "r": "tab:green", 

113 "i": "gold", 

114 "z": "tab:orange", 

115 "y": "tab:red", 

116 "N921": "tab:pink", 

117 }, 

118 ) 

119 vectorsToLoadList = ListField[str]( 

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

121 default=[ 

122 "visitId", 

123 "detector", 

124 "band", 

125 "ra", 

126 "decl", 

127 "zeroPoint", 

128 "psfSigma", 

129 "skyBg", 

130 "astromOffsetMean", 

131 "psfStarDeltaE1Median", 

132 "psfStarDeltaE2Median", 

133 "psfStarScaledDeltaSizeScatter", 

134 "psfTraceRadiusDelta", 

135 "llcra", 

136 "lrcra", 

137 "ulcra", 

138 "urcra", 

139 "llcdec", 

140 "lrcdec", 

141 "ulcdec", 

142 "urcdec", 

143 "xSize", 

144 "ySize", 

145 ], 

146 ) 

147 parametersToPlotList = ListField[str]( 

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

149 "plot these parameters for each band.", 

150 default=[ 

151 "psfSigma", 

152 "astromOffsetMean", 

153 "medianE", 

154 "psfStarScaledDeltaSizeScatter", 

155 "skyBg", 

156 "zeroPoint", 

157 ], 

158 ) 

159 tractsToPlotList = ListField[int]( 

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

161 default=None, 

162 optional=True, 

163 ) 

164 trimToTract = Field[bool]( 

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

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

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

168 default=False, 

169 ) 

170 doScatterInRaDec = Field[bool]( 

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

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

173 default=False, 

174 ) 

175 plotAllTractOutlines = Field[bool]( 

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

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

178 default=False, 

179 ) 

180 raDecLimitsDict = DictField[str, float]( 

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

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

183 "the keys raMin, ramax, decMin, decMax, e.g. " 

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

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

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

187 default=None, 

188 optional=True, 

189 ) 

190 plotDetectorOutline = Field[bool]( 

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

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

193 "in the inputs.", 

194 default=False, 

195 ) 

196 

197 def getInputSchema(self) -> KeyedDataSchema: 

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

199 for vector in self.vectorsToLoadList: 

200 base.append((vector, Vector)) 

201 return base 

202 

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

204 self._validateInput(data, **kwargs) 

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

206 

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

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

209 check that the data is consistent with Vector. 

210 """ 

211 needed = self.getFormattedInputSchema(**kwargs) 

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

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

214 }: 

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

216 for name, typ in needed: 

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

218 if isScalar and typ != Scalar: 

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

220 

221 if "medianE" in self.parametersToPlotList: 

222 if not all( 

223 vector in self.vectorsToLoadList 

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

225 ): 

226 raise RuntimeError( 

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

228 "psfStarDeltaE2Median must be included in vectorsToLoadList." 

229 ) 

230 if self.raDecLimitsDict is not None: 

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

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

233 raise RuntimeError( 

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

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

236 ) 

237 

238 def makePlot( 

239 self, 

240 data: KeyedData, 

241 plotInfo: Optional[Mapping[str, str]] = None, 

242 camera: Optional[Camera] = None, 

243 skymap: Optional[BaseSkyMap] = None, 

244 calibrateConfig: Optional[Config] = None, 

245 makeWarpConfig: Optional[Config] = None, 

246 **kwargs, 

247 ) -> Figure: 

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

249 

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

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

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

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

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

255 

256 Parameters 

257 ---------- 

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

259 The key-based catalog of data to plot. 

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

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

262 least) keys: 

263 

264 `"run"` 

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

266 `"tableName"` 

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

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

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

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

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

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

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

274 ``plotDetectorOutline`` is `True`). 

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

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

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

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

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

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

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

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

283 The persisted config used in the makeDirectWarp task 

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

285 in the run. 

286 

287 Returns 

288 ------- 

289 fig : `matplotlib.figure.Figure` 

290 The resulting figure. 

291 """ 

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

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

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

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

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

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

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

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

300 

301 if "medianE" in self.parametersToPlotList: 

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

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

304 ) 

305 

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

307 notFoundKeys = [] 

308 for zKey in self.parametersToPlotList: 

309 if zKey not in data.keys(): 

310 notFoundKeys.append(zKey) 

311 if len(notFoundKeys) > 0: 

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

313 

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

315 

316 if makeWarpConfig is None: 

317 maxEllipResidual = 0.007 

318 maxScaledSizeScatter = 0.019 

319 maxPsfTraceRadiusDelta = 0.7 

320 else: 

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

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

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

324 

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

326 if self.projection == "focalPlane": 

327 if camera is None: 

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

329 xKey = "xFp" 

330 yKey = "yFp" 

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

332 detFpDict = {} 

333 for det in camera: # type: ignore 

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

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

336 log.info("Number of SCIENCE detectors in {} camera = {}".format(cameraName, len(detFpDict))) 

337 xFpList = [] 

338 yFpList = [] 

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

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

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

342 

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

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

345 

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

347 xCorners, yCorners = zip(*corners) 

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

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

350 tractList: List[int] = [] 

351 elif self.projection == "raDec": 

352 xKey = "ra" 

353 yKey = "decl" 

354 xScatLen, yScatLen = 0, 0 

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

356 # RA/Dec entries. 

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

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

359 tractList = self.tractsToPlotList 

360 else: 

361 ras = data["ra"] 

362 decs = data["decl"] 

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

364 log.info("List of tracts overlapping data: {}".format(tractList)) 

365 tractLimitsDict = self._getTractLimitsDict(skymap, tractList) 

366 else: 

367 raise ValueError("Unknown projection: {}".format(self.projection)) 

368 

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

370 

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

372 data = self._trimDataToTracts(data, skymap) 

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

374 if nData == 0: 

375 raise RuntimeError( 

376 "No data to plot. Did your tract selection of " 

377 f"{self.tractsToPlotList} remove all data?" 

378 ) 

379 

380 if self.doScatterInRaDec: 

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

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

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

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

385 

386 dataDf = pd.DataFrame(data) 

387 nDataId = len(dataDf) 

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

389 if nDataId == 0: 

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

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

392 

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

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

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

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

397 if len(missingBandList) > 0: 

398 log.warning( 

399 "The band(s) {} are not included in self.sortedFullBandList. Please add them so " 

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

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

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

403 "(otherwise, if not present, it defaults to teal).".format(missingBandList) 

404 ) 

405 bandList.extend(missingBandList) 

406 log.info("Sorted list of existing bands: {}".format(bandList)) 

407 

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

409 if self.doScatterInRaDec: 

410 ptSize = min(3, 0.7 * ptSize) 

411 

412 nRow = len(bandList) 

413 nCol = len(self.parametersToPlotList) 

414 colMultiplier = 4 if nCol == 1 else 2.5 

415 fig, axes = plt.subplots( 

416 nRow, 

417 nCol, 

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

419 constrained_layout=True, 

420 ) 

421 

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

423 # all bands. 

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

425 vMinDict = {} 

426 vMaxDict = {} 

427 for zKey in self.parametersToPlotList: 

428 zKeySorted = dataDf[zKey].sort_values() 

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

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

431 if zKey == "medianE": 

432 vMaxDict[zKey] = maxEllipResidual 

433 elif zKey == "psfStarScaledDeltaSizeScatter": 

434 vMaxDict[zKey] = maxScaledSizeScatter 

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

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

437 elif zKey == "maxPsfTraceRadiusDelta": 

438 vMaxDict[zKey] = maxPsfTraceRadiusDelta 

439 else: 

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

441 

442 for iRow, band in enumerate(bandList): 

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

444 

445 nDataIdBand = len(dataBand) 

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

447 if nDataIdBand < 2: 

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

449 continue 

450 

451 for zKey in self.parametersToPlotList: 

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

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

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

455 zKeySorted = dataBand[zKey].sort_values() 

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

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

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

459 

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

461 if self.doScatterInRaDec: 

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

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

464 xScatter = scatRads * np.cos(scatTheta) 

465 yScatter = scatRads * np.sin(scatTheta) 

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

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

468 else: 

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

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

471 xLabel = xKey 

472 yLabel = yKey 

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

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

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

476 if self.projection == "focalPlane": 

477 for det in camera: # type: ignore 

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

479 detId = int(det.getId()) 

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

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

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

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

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

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

486 

487 if band not in self.bandLabelColorDict: 

488 log.warning( 

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

490 "to specify the label color (otherwise, it defaults to teal).".format(band) 

491 ) 

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

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

494 

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

496 if self.cmapName is None: 

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

498 else: 

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

500 

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

502 cmap.set_over("red") 

503 elif ( 

504 zKey in ["astromOffsetMean"] 

505 and self.projection != "raDec" 

506 and vMaxDict[zKey] >= maxMeanDistanceArcsec 

507 ): 

508 cmap.set_over("red") 

509 else: 

510 if self.cmapName is None: 

511 cmap.set_over("black") 

512 else: 

513 cmap.set_over("lemonchiffon") 

514 

515 if zKey in ["psfSigma"]: 

516 cmap.set_under("red") 

517 else: 

518 if self.cmapName is None: 

519 cmap.set_under("lemonchiffon") 

520 else: 

521 cmap.set_under("black") 

522 

523 titleStr = "band: {} nVisit: {} nData: {}".format(band, nVisitBand, nDataIdBand) 

524 

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

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

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

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

529 ax.set_aspect("equal") 

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

531 

532 if self.projection == "focalPlane": 

533 for det in camera: # type: ignore 

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

535 corners = det.getCorners(FOCAL_PLANE) 

536 xCorners, yCorners = zip(*corners) 

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

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

539 detId = int(det.getId()) 

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

541 if len(perDetData) < 1: 

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

543 continue 

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

545 log.debug( 

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

547 ) 

548 continue 

549 pcm = plotProjectionWithBinning( 

550 ax, 

551 perDetData["xScat"].values, 

552 perDetData["yScat"].values, 

553 perDetData[zKey].values, 

554 cmap, 

555 xFpMin, 

556 xFpMax, 

557 yFpMin, 

558 yFpMax, 

559 xNumBins=1, 

560 fixAroundZero=False, 

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

562 isSorted=False, 

563 vmin=vMinDict[zKey], 

564 vmax=vMaxDict[zKey], 

565 showExtremeOutliers=self.showExtremeOutliers, 

566 scatPtSize=ptSize, 

567 ) 

568 if self.projection == "raDec": 

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

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

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

572 pcm = plotProjectionWithBinning( 

573 ax, 

574 dataBand["xScat"].values, 

575 dataBand["yScat"].values, 

576 dataBand[zKey].values, 

577 cmap, 

578 raMin, 

579 raMax, 

580 decMin, 

581 decMax, 

582 xNumBins=self.nBins, 

583 fixAroundZero=False, 

584 nPointBinThresh=self.nPointBinThresh, 

585 isSorted=False, 

586 vmin=vMinDict[zKey], 

587 vmax=vMaxDict[zKey], 

588 showExtremeOutliers=self.showExtremeOutliers, 

589 scatPtSize=ptSize, 

590 ) 

591 

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

593 # equal axis ratio panels. 

594 if iRow == 0 and iCol == 0: 

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

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

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

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

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

600 else: 

601 xLimMin, xLimMax = ax.get_xlim() 

602 yLimMin, yLimMax = ax.get_ylim() 

603 if self.projection == "focalPlane": 

604 minDim = ( 

605 max( 

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

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

608 ) 

609 / 2 

610 ) 

611 xLimMin = min(-minDim, xLimMin) 

612 xLimMax = max(minDim, xLimMax) 

613 if self.trimToTract: 

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

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

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

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

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

619 xDelta = xLimMax - xLimMin 

620 xLimMin -= 0.04 * xDelta 

621 xLimMax += 0.04 * xDelta 

622 

623 xLimMin, xLimMax, yLimMin, yLimMax = _setLimitsToEqualRatio( 

624 xLimMin, xLimMax, yLimMin, yLimMax 

625 ) 

626 limRange = xLimMax - xLimMin 

627 yTickFmt = _tickFormatter(yLimMin * 10) 

628 

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

630 allTractsList = [] 

631 for tractInfo in skymap: # type: ignore 

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

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

634 if ( 

635 centerRa > xLimMin 

636 and centerRa < xLimMax 

637 and centerDec > yLimMin 

638 and centerDec < yLimMax 

639 ): 

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

641 tractLimitsDict = self._getTractLimitsDict(skymap, allTractsList) 

642 

643 upperHandles = [] 

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

645 patch = mpl.patches.Circle( 

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

647 radius=scatRad, 

648 facecolor="gray", 

649 edgecolor="None", 

650 alpha=0.2, 

651 label="scatter area", 

652 ) 

653 ax.add_patch(patch) 

654 upperHandles.append(patch) 

655 

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

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

658 if camera is None: 

659 if iRow == 0 and iCol == 0: 

660 log.warning( 

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

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

663 ) 

664 else: 

665 # Calculate area of polygon with known vertices. 

666 x1, x2, x3, x4 = ( 

667 dataBand["llcra"], 

668 dataBand["lrcra"], 

669 dataBand["urcra"], 

670 dataBand["ulcra"], 

671 ) 

672 y1, y2, y3, y4 = ( 

673 dataBand["llcdec"], 

674 dataBand["lrcdec"], 

675 dataBand["urcdec"], 

676 dataBand["ulcdec"], 

677 ) 

678 areaDeg = ( 

679 np.abs( 

680 (x1 * y2 - y1 * x2) 

681 + (x2 * y3 - y2 * x3) 

682 + (x3 * y4 - y3 * x4) 

683 + (x4 * y1 - y4 * x1) 

684 ) 

685 / 2.0 

686 ) 

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

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

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

690 

691 patch = mpl.patches.Rectangle( 

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

693 detWidthDeg, 

694 detHeightDeg, 

695 facecolor="turquoise", 

696 alpha=0.3, 

697 label="det size", 

698 ) 

699 ax.add_patch(patch) 

700 upperHandles.append(patch) 

701 

702 if self.projection == "raDec": 

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

704 ax.legend( 

705 handles=upperHandles, 

706 loc="upper left", 

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

708 edgecolor="black", 

709 framealpha=0.4, 

710 fontsize=5, 

711 ) 

712 

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

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

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

716 if len(tractList) > 0: 

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

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

719 if limRange <= 5 * tractWidthDeg: 

720 deltaLim = 0.05 * limRange 

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

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

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

724 if ( 

725 centerRa > xLimMin + deltaLim 

726 and centerRa < xLimMax - deltaLim 

727 and centerDec > yLimMin + deltaLim 

728 and centerDec < yLimMax - deltaLim 

729 ): 

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

731 fontSize = 3 if limRange < 20 else 2 

732 ax.annotate( 

733 str(tract), 

734 tractLimits["center"], 

735 va="center", 

736 ha="center", 

737 color="dimgray", 

738 fontsize=fontSize, 

739 annotation_clip=True, 

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

741 ) 

742 if self.projection == "raDec": 

743 ax.set_xlim(xLimMax, xLimMin) 

744 else: 

745 ax.set_xlim(xLimMin, xLimMax) 

746 ax.set_ylim(yLimMin, yLimMax) 

747 ax.yaxis.set_major_formatter(yTickFmt) 

748 

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

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

751 # has the same effect on all colorbars. 

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