Coverage for python / lsst / analysis / tools / actions / plot / propertyMapPlot.py: 9%

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__ = ( 

24 "PerTractPropertyMapPlot", 

25 "SurveyWidePropertyMapPlot", 

26) 

27 

28import logging 

29from typing import Mapping, Union 

30 

31import lsst.pex.config as pexConfig 

32import lsst.sphgeom as sphgeom 

33import matplotlib.patheffects as mpl_path_effects 

34import numpy as np 

35import skyproj 

36from healsparse.healSparseMap import HealSparseMap 

37from lsst.analysis.tools.tasks.propertyMapAnalysis import ( 

38 PerTractPropertyMapAnalysisConfig, 

39 SurveyWidePropertyMapAnalysisConfig, 

40) 

41from lsst.skymap.tractInfo import ExplicitTractInfo 

42from lsst.utils.plotting import make_figure, set_rubin_plotstyle 

43from matplotlib import cm, rc_context 

44from matplotlib.figure import Figure 

45from matplotlib.legend_handler import HandlerTuple 

46from mpl_toolkits.axes_grid1 import make_axes_locatable 

47from mpl_toolkits.axes_grid1.axes_size import AxesX, AxesY, Fraction 

48 

49from ...interfaces import KeyedData, PlotAction 

50 

51_LOG = logging.getLogger(__name__) 

52 

53# Holds unit renames to match style guidelines: {"old_unit": "new_unit"}. 

54unitRenameDict = {"mag(AB)": r"$\rm mag_{AB}$"} 

55 

56 

57def getZoomedExtent(box, n): 

58 """Get a zoomed plot extent in degrees for a bounding box. 

59 

60 Parameters 

61 ---------- 

62 box : `lsst.sphgeom.NormalizedAngleInterval` 

63 The bounding box to get an extent for. 

64 n : `float` 

65 Zoom factor; for instance, n=2 means zooming in 2 times at the 

66 center. Must be positive. 

67 

68 Returns 

69 ------- 

70 `tuple` [`float`] 

71 New extent as (new_lon_min, new_lon_max, new_lat_min, new_lat_max). 

72 """ 

73 lon = box.getLon() 

74 lat = box.getLat() 

75 if n > 1: 

76 factor = (1.0 - 1.0 / n) / 2.0 

77 lon.erodeBy(sphgeom.NormalizedAngle(factor * box.getWidth())) 

78 lat.erodeBy(sphgeom.NormalizedAngle(factor * box.getHeight())) 

79 elif n < 1: 

80 factor = (1.0 / n - 1) / 2.0 

81 lon.dilateBy(sphgeom.NormalizedAngle(factor * box.getWidth())) 

82 lat.dilateBy(sphgeom.NormalizedAngle(factor * box.getHeight())) 

83 

84 extent = ( 

85 lon.getA().asDegrees(), 

86 lon.getB().asDegrees(), 

87 lat.getA().asDegrees(), 

88 lat.getB().asDegrees(), 

89 ) 

90 return extent 

91 

92 

93def getLongestSuffixMatch(s, options): 

94 """Find the longest suffix in the provided list that matches the end of 

95 the given string. 

96 

97 Parameters 

98 ---------- 

99 s : `str` 

100 The target string for which we want to find a matching suffix. 

101 options : `list` [`str`] 

102 A list of potential suffix strings to match against the target 

103 string `s`. 

104 

105 Returns 

106 ------- 

107 `str` 

108 The longest matching suffix from the `options` list. If no match is 

109 found, returns `None`. 

110 """ 

111 return next((opt for opt in sorted(options, key=len, reverse=True) if s.endswith(opt)), None) 

112 

113 

114def addTextToColorbar( 

115 cb, text, orientation="vertical", color="black", fontsize=14, fontweight="bold", alpha=0.8 

116): 

117 """Helper method to add text inside the horizontal colorbar. 

118 

119 Parameters 

120 ---------- 

121 cb : `~matplotlib.colorbar.Colorbar` 

122 The colorbar object. 

123 text : `str` 

124 The text to add. 

125 orientation : `str`, optional 

126 The orientation of the colorbar. Can be either "vertical" or 

127 "horizontal". 

128 fontsize : `int`, optional 

129 The fontsize of the text. 

130 fontweight : `str`, optional 

131 The fontweight of the text. 

132 alpha : `float`, optional 

133 The alpha value of the text. 

134 

135 Returns 

136 ------- 

137 `None` 

138 The text is added to the colorbar in place. 

139 """ 

140 if color is None: 

141 color = "black" 

142 vmid = (cb.vmin + cb.vmax) / 2 

143 positions = {"vertical": (0.5, vmid), "horizontal": (vmid, 0.5)} 

144 cbtext = cb.ax.text( 

145 *positions[orientation], 

146 text, 

147 color=color, 

148 va="center", 

149 ha="center", 

150 fontsize=fontsize, 

151 fontweight=fontweight, 

152 rotation=orientation, 

153 alpha=alpha, 

154 ) 

155 # Add a distinct outline around the text for better visibility in 

156 # various backgrounds. 

157 cbtext.set_path_effects( 

158 [mpl_path_effects.Stroke(linewidth=4, foreground="white", alpha=0.8), mpl_path_effects.Normal()] 

159 ) 

160 

161 

162class CustomHandler(HandlerTuple): 

163 """Custom legend handler to overlay multiple patches for a single 

164 legend entry. 

165 

166 This handler class inherits from `HandlerTuple` and is designed to 

167 handle cases where multiple artists (e.g., patches) need to be overlaid 

168 on top of each other in a single legend entry, as opposed to 

169 side-by-side which is the default behavior of `HandlerTuple`. 

170 

171 Methods 

172 ------- 

173 create_artists: 

174 Override the `create_artists` method of `HandlerTuple` to modify 

175 the positioning of the artists so that they overlay directly on top 

176 of one another in the legend. 

177 

178 Example 

179 ------- 

180 # Plot some data. 

181 line, = ax.plot(x, y, label="Sample Line") 

182 

183 # Use CustomHandler for overlaid patches and also include the regular 

184 # line legend if desired. 

185 handles = [(patch1, patch2), line] 

186 labels = ['Overlaid Patches', line.get_label()] 

187 leg = ax.legend( 

188 handles, labels, handler_map={tuple: CustomHandler()}, loc="best" 

189 ) 

190 """ 

191 

192 def create_artists(self, legend, orig_handle, xdescent, ydescent, width, height, fontsize, trans): 

193 artists = HandlerTuple.create_artists( 

194 self, legend, orig_handle, xdescent, ydescent, width, height, fontsize, trans 

195 ) 

196 # Overlay the two patches. 

197 for a in artists: 

198 a.set_transform(trans) 

199 return artists 

200 

201 

202class PerTractPropertyMapPlot(PlotAction): 

203 draw_patch_bounds = pexConfig.Field[bool]( 

204 doc="Whether to draw patch inner boundaries or not", 

205 default=False, 

206 ) 

207 label_patches = pexConfig.ChoiceField[str]( 

208 doc="Which patches to label by ID", 

209 default="none", 

210 allowed={ 

211 "all": "All patches", 

212 "edge": "Edge patches only", 

213 "none": "No labels", 

214 }, 

215 ) 

216 plotName = pexConfig.Field[str](doc="The name for the plotting task.", optional=True) 

217 

218 def __call__( 

219 self, 

220 data: KeyedData, 

221 tractInfo: ExplicitTractInfo, 

222 plotConfig: PerTractPropertyMapAnalysisConfig, 

223 plotInfo: Mapping[str, Union[Mapping[str, str], str, int]], 

224 **kwargs, 

225 ) -> Mapping[str, Figure]: 

226 self._validateInput(data, tractInfo, plotConfig, plotInfo) 

227 return self.makePlot(data, tractInfo, plotConfig, plotInfo) 

228 

229 def _validateInput( 

230 self, 

231 data: KeyedData, 

232 tractInfo: ExplicitTractInfo, 

233 plotConfig: PerTractPropertyMapAnalysisConfig, 

234 plotInfo: Mapping[str, Union[Mapping[str, str], str, int]], 

235 ) -> None: 

236 """Validate the input data.""" 

237 

238 if not isinstance(tractInfo, ExplicitTractInfo): 

239 raise TypeError(f"Input `tractInfo` type must be {ExplicitTractInfo} not {type(tractInfo)}.") 

240 

241 if not isinstance(plotConfig, PerTractPropertyMapAnalysisConfig): 

242 raise TypeError( 

243 "`plotConfig` must be a " 

244 "`lsst.analysis.tools.tasks.propertyMapAnalysis.PerTractPropertyMapAnalysisConfig`. " 

245 f"Got {type(plotConfig)}." 

246 ) 

247 

248 if not isinstance(plotInfo, dict): 

249 raise TypeError("`plotConfig` must be a dictionary.") 

250 

251 if not plotConfig.publicationStyle: 

252 zoomFactors = plotConfig.zoomFactors 

253 isListOfFloats = isinstance(zoomFactors, pexConfig.listField.List) and all( 

254 isinstance(zf, float) for zf in zoomFactors 

255 ) 

256 if not (isListOfFloats and len(zoomFactors) == 2) or any(zf <= 1 for zf in zoomFactors): 

257 raise TypeError( 

258 "`zoomFactors` must be a two-element `lsst.pex.config.listField.List` of floats > 1." 

259 ) 

260 

261 for atool in plotConfig.atools: 

262 if not isinstance(atool.nBinsHist, int) or atool.nBinsHist <= 0: 

263 raise ValueError( 

264 f"`nBinsHist` for property `{atool.process.buildActions.data.mapKey}` must be a positive " 

265 f"integer. Got {atool.nBinsHist}." 

266 ) 

267 

268 # Identify any invalid entries in `data`. 

269 invalidEntries = { 

270 key: pytype 

271 for key, pytype in {k: v.ref.datasetType.storageClass.pytype for k, v in data.items()}.items() 

272 if pytype != HealSparseMap 

273 } 

274 

275 # If any invalid entries are found, raise a TypeError with details. 

276 if invalidEntries: 

277 errorMessage = "; ".join( 

278 f"`{key}` should be {HealSparseMap}, got {type_}" for key, type_ in invalidEntries.items() 

279 ) 

280 raise TypeError(f"Invalid input types found in `data`: {errorMessage}") 

281 

282 def addPlotInfo( 

283 self, 

284 fig: Figure, 

285 plotInfo: Mapping[str, Union[Mapping[str, str], str, int]], 

286 toolName: str, 

287 ) -> Figure: 

288 """Add useful information to the plot. 

289 

290 Parameters 

291 ---------- 

292 fig : `matplotlib.figure.Figure` 

293 The figure to add the information to. 

294 plotInfo : `dict` 

295 A dictionary of the plot information. 

296 toolName : `str` 

297 The name of the tool used to generate the plot. 

298 

299 Returns 

300 ------- 

301 fig : `matplotlib.figure.Figure` 

302 The figure with the information added. 

303 """ 

304 

305 run = plotInfo["run"] 

306 tableType = f"\nTable: {plotInfo['tableNames'][toolName]}" 

307 

308 dataIdText = f"Tract: {plotInfo['tract']}, Band: {plotInfo['band']}" 

309 propertyDescription = plotInfo["description"] 

310 # Lowercase the first letter unless the string starts with more than 

311 # one uppercase letter in which case we keep it as is, e.g. PSF, DCR. 

312 if propertyDescription[0].isupper() and propertyDescription[1].islower(): 

313 propertyDescription = propertyDescription[0].lower() + propertyDescription[1:] 

314 mapText = ( 

315 f", Property: {propertyDescription}, " 

316 f"Unit: {plotInfo['unit']}, " 

317 f"Operation: {plotInfo['operation']}, " 

318 f"Coadd: {plotInfo['coaddName']}" 

319 ) 

320 geomText = ( 

321 f", Valid area: {plotInfo['valid_area']:.2f} sq. deg., " 

322 + f"NSIDE: {plotInfo['nside']}, projection: {plotInfo['projection']}" 

323 ) 

324 infoText = f"\n{dataIdText}{mapText}" 

325 

326 fig.text( 

327 0.04, 

328 0.965, 

329 f'{plotInfo["plotName"]} of {plotInfo["property"]}', 

330 fontsize=19, 

331 transform=fig.transFigure, 

332 ha="left", 

333 va="top", 

334 ) 

335 t = fig.text( 

336 0.04, 

337 0.942, 

338 f"{run}{tableType}{geomText}{infoText}", 

339 fontsize=15, 

340 transform=fig.transFigure, 

341 alpha=0.6, 

342 ha="left", 

343 va="top", 

344 ) 

345 t.set_linespacing(1.4) 

346 

347 return fig 

348 

349 def makePlot( 

350 self, 

351 data: KeyedData, 

352 tractInfo: ExplicitTractInfo, 

353 plotConfig: PerTractPropertyMapAnalysisConfig, 

354 plotInfo: Mapping[str, Union[Mapping[str, str], str, int]], 

355 ) -> Mapping[str, Figure]: 

356 """Make the survey property map plot. 

357 

358 Parameters 

359 ---------- 

360 data : `KeyedData` 

361 The HealSparseMap to plot the points from. 

362 tractInfo: `~lsst.skymap.tractInfo.ExplicitTractInfo` 

363 The tract info object. 

364 plotConfig : 

365 `~lsst.analysis.tools.tasks.perTractPropertyMapAnalysis. 

366 PerTractPropertyMapAnalysisConfig` 

367 The configuration for the plot. 

368 plotInfo : `dict` 

369 A dictionary of information about the data being plotted. 

370 

371 Returns 

372 ------- 

373 fig : `~matplotlib.figure.Figure` 

374 The resulting figure. 

375 """ 

376 

377 set_rubin_plotstyle() 

378 

379 # 'plotName' by default is constructed from the attribute specified in 

380 # 'atools.<attribute>' in the pipeline YAML. If it is explicitly 

381 # set in `~lsst.analysis.tools.atools.propertyMap.PropertyMapTool`, 

382 # it will override this default. 

383 if self.plotName: 

384 # Set the plot name using 'produce.plot.plotName' from 

385 # PropertyMapTool's instance. 

386 plotInfo["plotName"] = self.plotName 

387 

388 # Plotting customization. 

389 colorbarLocation = "right" 

390 colorbarTickLabelSize = 14 

391 colorBarAspect = 16 

392 colorBarFraction = 0.07 # w.r.t. the axes size. 

393 colorBarPad = 0 # w.r.t. the axes size. 

394 rcparams = { 

395 "axes.labelsize": 22 if plotConfig.publicationStyle else 19, 

396 "axes.linewidth": 1.8, 

397 "xtick.labelsize": 15 if plotConfig.publicationStyle else 14, 

398 "ytick.labelsize": 15 if plotConfig.publicationStyle else 14, 

399 } 

400 

401 colorbarKwargs = dict(plotConfig.colorbarKwargs) 

402 projectionKwargs = dict(plotConfig.projectionKwargs) 

403 

404 if plotConfig.publicationStyle: 

405 zoomFactors = [] 

406 _LOG.info( 

407 "Zoom factors are not used in publication-style plots. " 

408 "Only the full-tract map will be plotted." 

409 ) 

410 if "cmap" in colorbarKwargs and colorbarKwargs["cmap"] != "viridis": 

411 _LOG.warning( 

412 "Color map set to 'viridis' for publication style plots. " 

413 f"The color map '{colorbarKwargs['cmap']}' set in the config will be ignored." 

414 ) 

415 colorbarKwargs["cmap"] = "viridis" 

416 cmap = cm.get_cmap(colorbarKwargs["cmap"]) 

417 # Colorbar text color only, not used for any histogram. 

418 histColors = [cmap(0.1)] 

419 labelpad = 14 

420 else: 

421 zoomFactors = plotConfig.zoomFactors 

422 colorbarKwargs["cmap"] = colorbarKwargs.get("cmap", "viridis") 

423 # Muted green for the full map, and muted red and blue for the two 

424 # zoomed-in maps. Used for boxes, colorbar texts and histograms. 

425 histColors = ["#265D40", "#8B0000", "#00008B"] 

426 labelpad = 11 

427 

428 toolName = data["data"].ref.datasetType.name 

429 mapName = toolName.replace("_map_", "_") 

430 mapData = data["data"].get() 

431 

432 with rc_context(rcparams): 

433 if plotConfig.publicationStyle: 

434 # Make a single plot for the full tract. 

435 fig = make_figure(figsize=(8, 8)) 

436 ax1 = fig.add_subplot(111) 

437 fig.subplots_adjust(left=0.133, right=0.888, bottom=0.1, top=0.855) 

438 else: 

439 # Make a 2x2 grid of subplots for the full tract, two zoomed-in 

440 # views, and a histogram of values. 

441 fig = make_figure(figsize=(16, 16)) 

442 ax1 = fig.add_subplot(221) 

443 ax2 = fig.add_subplot(222) 

444 ax3 = fig.add_subplot(223) 

445 ax4 = fig.add_subplot(224) 

446 

447 # Reduce whitespace but leave some room at the top for info. 

448 fig.subplots_adjust(left=0.064, right=0.96, top=0.855, bottom=0.07, wspace=0.275, hspace=0.24) 

449 

450 # Get the values for the valid pixels of the full tract. 

451 values = mapData[mapData.valid_pixels] 

452 goodValues = np.isfinite(values) 

453 values = values[goodValues] # As a precaution. 

454 

455 # Make a concise human-readable label for the plot. 

456 plotInfo["unit"] = "N/A" # Unless overridden. 

457 if hasattr(mapData, "metadata") and all( 

458 key in mapData.metadata for key in ["DESCRIPTION", "OPERATION", "UNIT"] 

459 ): 

460 hasMetadata = True 

461 metadata = mapData.metadata 

462 plotInfo["description"] = metadata["DESCRIPTION"] 

463 plotInfo["operation"] = metadata["OPERATION"] 

464 if metadata["UNIT"]: 

465 plotInfo["unit"] = metadata["UNIT"] 

466 elif metadata["UNIT"] == "": 

467 plotInfo["unit"] = "dimensionless" 

468 else: 

469 hasMetadata = False 

470 plotInfo["operation"] = getLongestSuffixMatch( 

471 mapName, ["min", "max", "mean", "weighted_mean", "sum"] 

472 ).replace("_", " ") 

473 plotInfo["coaddName"] = mapName.split("Coadd_")[0] 

474 plotInfo["operation"] = plotInfo["operation"].replace("minimum", "min").replace("maximum", "max") 

475 propertyName = mapName[len(f"{plotInfo['coaddName']}Coadd_") : -len(plotInfo["operation"])].strip( 

476 "_" 

477 ) 

478 if not hasMetadata: 

479 # Infer the property description from the map name (all 

480 # lower case), and properly handle formatting. 

481 plotInfo["description"] = ( 

482 propertyName.replace("_", " ") 

483 .replace("psf", "PSF") 

484 .replace("dcr", "DCR") 

485 .replace("dra", r"$\Delta$RA") 

486 .replace("ddec", r"$\Delta$Dec") 

487 ) 

488 plotInfo["property"] = ( 

489 propertyName.replace("_", " ") 

490 .title() # Capitalize and handle edge cases below. 

491 .replace("Psf", "PSF") 

492 .replace("Dcr", "DCR") 

493 .replace("Dra", r"$\Delta$RA") 

494 .replace("Ddec", r"$\Delta$Dec") 

495 .replace("E1", "e1") 

496 .replace("E2", "e2") 

497 ) 

498 

499 # Handle unit renaming. 

500 if plotInfo["unit"] in unitRenameDict: 

501 plotInfo["unit"] = unitRenameDict[plotInfo["unit"]] 

502 

503 atool = getattr(plotConfig.atools, toolName) 

504 nBinsHist = atool.nBinsHist 

505 

506 # Use the full tract bounding box to set the default extent. 

507 box = tractInfo.getOuterSkyPolygon().getBoundingBox() 

508 ctr_lonlat = box.getCenter() 

509 ctr_lon = ctr_lonlat.getLon().asDegrees() 

510 ctr_lat = ctr_lonlat.getLat().asDegrees() 

511 

512 # Prepare plot elements for the full tract and optional zoomed-in 

513 # views. 

514 if plotConfig.publicationStyle: 

515 skyprojAxes = [ax1] 

516 else: 

517 skyprojAxes = [ax1, ax3, ax4] 

518 

519 zoomIdx = [] 

520 for ax, zoomFactor, histColor in zip(skyprojAxes, [1.0, *zoomFactors], histColors): 

521 extent = getZoomedExtent(box, zoomFactor) 

522 

523 sp = skyproj.GnomonicSkyproj( 

524 ax=ax, 

525 lon_0=ctr_lon, 

526 lat_0=ctr_lat, 

527 extent=extent, 

528 **projectionKwargs, 

529 ) 

530 

531 sp.draw_hspmap(mapData, zoom=False, cmap=colorbarKwargs["cmap"]) 

532 

533 sp.ax.set_xlabel("R.A.", labelpad=labelpad, fontsize=rcparams["axes.labelsize"]) 

534 sp.ax.set_ylabel("Dec.", labelpad=labelpad, fontsize=rcparams["axes.labelsize"]) 

535 

536 if self.draw_patch_bounds: 

537 label_all = self.label_patches == "all" 

538 label_edge = self.label_patches == "edge" 

539 if label_all or label_edge: 

540 patchids_max = tuple(x - 1 for x in tractInfo.getNumPatches()) 

541 lon_min, lon_max = min(extent[:2]), max(extent[:2]) 

542 lat_min, lat_max = min(extent[2:]), max(extent[2:]) 

543 

544 for patchInfo in tractInfo: 

545 vertices = patchInfo.getInnerSkyPolygon().getVertices() 

546 clipped = tractInfo.inner_sky_region.clipTo( 

547 sphgeom.Box(sphgeom.LonLat(vertices[0]), sphgeom.LonLat(vertices[2])) 

548 ) 

549 lonlats_patch = np.array( 

550 [ 

551 [x.asDegrees() for x in (lonlat.getA(), lonlat.getB())] 

552 for lonlat in (clipped.getLon(), clipped.getLat()) 

553 ] 

554 ) 

555 lons_patch = np.concat((lonlats_patch[0, :], lonlats_patch[0, ::-1])) 

556 lats_patch = np.repeat(lonlats_patch[1, :], 2) 

557 

558 sp.draw_polygon(lons_patch, lats_patch, edgecolor="gray") 

559 label_id = label_all or ( 

560 label_edge 

561 and ( 

562 (patchInfo.index[0] == 0) 

563 or (patchInfo.index[1] == 0) 

564 or (patchInfo.index[0] == patchids_max[0]) 

565 or (patchInfo.index[1] == patchids_max[1]) 

566 ) 

567 ) 

568 if label_id: 

569 lon_label = clipped.getCenter().getLon().asDegrees() 

570 lat_label = clipped.getCenter().getLat().asDegrees() 

571 if (lon_min < lon_label < lon_max) and (lat_min < lat_label < lat_max): 

572 sp.ax.text( 

573 lon_label, 

574 lat_label, 

575 patchInfo.getSequentialIndex(), 

576 ha="center", 

577 va="center", 

578 size=7, 

579 c=[0, 0, 0, 0.5], 

580 ) 

581 

582 # Specify the size and padding of the colorbar axes with 

583 # respect to the main axes. 

584 refAx = AxesX(ax) if colorbarLocation in ("left", "right") else AxesY(ax) 

585 cbsize = Fraction(colorBarFraction, refAx) 

586 cbpad = Fraction(colorBarPad, refAx) 

587 

588 # Make divider and a colorbar axes to be attached to the main 

589 # axes. 

590 divider = make_axes_locatable(ax) 

591 cax = divider.append_axes(colorbarLocation, size=cbsize, pad=cbpad) 

592 

593 cbar = sp.draw_colorbar( 

594 **{ 

595 "location": colorbarLocation, 

596 "aspect": colorBarAspect, 

597 "fraction": colorBarFraction, 

598 "pad": colorBarPad, 

599 "cax": cax, 

600 **colorbarKwargs, 

601 } 

602 ) 

603 cbar.ax.tick_params(labelsize=colorbarTickLabelSize) 

604 if plotConfig.publicationStyle: 

605 if plotInfo["unit"] not in ["dimensionless", "N/A"]: 

606 cbarText = f"{plotInfo['property']} ({plotInfo['unit']})" 

607 else: 

608 cbarText = plotInfo["property"] 

609 else: 

610 cbarText = ( 

611 "Full Tract" if zoomFactor == 1.0 else f"{self.prettyPrintFloat(zoomFactor)}x Zoom" 

612 ) 

613 addTextToColorbar(cbar, cbarText, color=histColor) 

614 if zoomFactor == 1.0: 

615 # Store the "full tract" map so that we can overplot 

616 # the zoom rectangles. 

617 spf = sp 

618 else: 

619 # Create a rectangle for the zoomed-in region. 

620 x0, x1, y0, y1 = extent 

621 for c, ls, lw in zip(["white", histColor], ["solid", "dashed"], [3.5, 1.5]): 

622 spf.draw_polygon( 

623 [x0, x0, x1, x1], 

624 [y0, y1, y1, y0], 

625 facecolor="none", 

626 edgecolor=c, 

627 linestyle=ls, 

628 linewidth=lw, 

629 alpha=0.8, 

630 ) 

631 zoomText = spf.ax.text( 

632 (x0 + x1) / 2, 

633 y0, 

634 f"{self.prettyPrintFloat(zoomFactor)}x", 

635 color=histColor, 

636 fontsize=14, 

637 fontweight="bold", 

638 alpha=0.8, 

639 ha="center", 

640 va="bottom", 

641 ) 

642 # Add a distinct outline around the text for better 

643 # visibility in various backgrounds. 

644 zoomText.set_path_effects( 

645 [ 

646 mpl_path_effects.Stroke(linewidth=4, foreground="white", alpha=0.8), 

647 mpl_path_effects.Normal(), 

648 ] 

649 ) 

650 # Get the indices of pixels in the zoomed-in region. 

651 pos = mapData.valid_pixels_pos() 

652 # Reversed axes consideration. 

653 xmin, xmax = sorted([x0, x1]) 

654 idx = (pos[0] > xmin) & (pos[0] < xmax) & (pos[1] > y0) & (pos[1] < y1) 

655 zoomIdx.append(idx[goodValues]) 

656 

657 # Calculate weights for each bin to ensure that the peak of the 

658 # histogram reaches 1. 

659 if len(values) == 0 or np.ptp(values) < 1e-12: 

660 # np.histogram cannot use 100 bins if no variance in values 

661 weights = np.ones_like(values) / values.size 

662 nBinsHist = 1 

663 _LOG.info(f"No variance in {toolName}; set histogram to a single bin.") 

664 else: 

665 weights = np.ones_like(values) / np.histogram(values, bins=nBinsHist)[0].max() 

666 

667 if not plotConfig.publicationStyle: 

668 # Compute full-tract histogram and get its bins. 

669 # NOTE: `exposure_time` histograms are quantized and look more 

670 # bar-like, so they look better with fewer bins. 

671 bins = ax2.hist( 

672 values, 

673 bins=nBinsHist, 

674 label="Full Tract", 

675 color=histColors[0], 

676 weights=weights, 

677 alpha=0.7, 

678 )[1] 

679 

680 # Align the histogram (top right panel) with the skyproj plots. 

681 pos1 = spf.ax.get_position() # Top left. 

682 pos4 = sp.ax.get_position() # Bottom right. 

683 cbarWidth = cbar.ax.get_position().height / colorBarAspect 

684 # NOTE: cbarWidth != cbar.ax.get_position().width 

685 ax2.set_position([pos4.x0, pos1.y0, pos4.width + cbarWidth, pos1.height]) 

686 

687 # Overplot the histograms for the zoomed-in plots. 

688 for zoomFactor, zidx, color, linestyle, hatch in zip( 

689 zoomFactors, zoomIdx, histColors[1:], ["solid", "dotted"], ["//", "xxxx"] 

690 ): 

691 weights = np.ones_like(values[zidx]) / np.histogram(values[zidx], bins=bins)[0].max() 

692 histLabel = f"{self.prettyPrintFloat(zoomFactor)}x Zoom" 

693 histValues = ax2.hist( 

694 values[zidx], 

695 bins=bins, 

696 label=histLabel, 

697 color=color, 

698 weights=weights, 

699 histtype="step", 

700 linewidth=2, 

701 linestyle=linestyle, 

702 alpha=0.6, 

703 )[0] 

704 # Fill the area under the step. 

705 ax2.fill_between( 

706 bins[:-1], 

707 histValues, 

708 step="post", 

709 color=color, 

710 alpha=0.2, 

711 hatch=hatch, 

712 label="hidden", 

713 ) 

714 

715 # Set labels and legend. 

716 xlabel = plotInfo["property"] 

717 if plotInfo["unit"] not in ["dimensionless", "N/A"]: 

718 xlabel += f" ({plotInfo['unit']})" 

719 xtext = ax2.set_xlabel(xlabel, labelpad=labelpad) 

720 ytext = ax2.set_ylabel("Normalized Count", labelpad=labelpad) 

721 xtext.set_fontsize(rcparams["axes.labelsize"]) 

722 ytext.set_fontsize(rcparams["axes.labelsize"]) 

723 

724 # Get handles and labels from the axis. 

725 handles, labels = ax2.get_legend_handles_labels() 

726 

727 # Add a legend with custom handler that combines the handle 

728 # pairs for the zoomed-in cases. 

729 handles = [handles[0], (handles[1], handles[2]), (handles[3], handles[4])] 

730 while "hidden" in labels: 

731 labels.remove("hidden") 

732 legend = ax2.legend( 

733 handles, 

734 labels, 

735 handler_map={tuple: CustomHandler()}, 

736 loc="best", 

737 frameon=False, 

738 fontsize=15, 

739 ) 

740 

741 for line, text in zip(handles, legend.get_texts()): 

742 if isinstance(line, tuple): 

743 # Use the first handle to get the color. 

744 line = line[0] 

745 color = line.get_edgecolor() if line.get_facecolor()[-1] == 0 else line.get_facecolor() 

746 text.set_color(color) 

747 

748 # Add extra info to plotInfo. 

749 plotInfo["projection"] = "Gnomonic" 

750 plotInfo["nside"] = mapData.nside_sparse 

751 plotInfo["valid_area"] = mapData.get_valid_area() 

752 

753 # Add useful information to the plot. 

754 self.addPlotInfo(fig, plotInfo, toolName) 

755 style = "" 

756 else: 

757 style = "publication-style " 

758 fig.suptitle( 

759 f"{plotInfo['description']} {plotInfo['operation']} map", 

760 fontsize=18.5, 

761 ha="center", 

762 va="top", 

763 y=0.985, 

764 ) 

765 fig.text( 

766 0.5, 

767 0.925, 

768 f"Tract: {plotInfo['tract']}, Band: {plotInfo['band']}, Coadd: {plotInfo['coaddName']}", 

769 ha="center", 

770 fontsize=15.5, 

771 ) 

772 

773 _LOG.info( 

774 f"Made {style}per-tract property map plot for dataset type '{toolName}', " 

775 f"tract: {plotInfo['tract']}, band: '{plotInfo['band']}'." 

776 ) 

777 

778 return fig 

779 

780 @staticmethod 

781 def prettyPrintFloat(n): 

782 if n.is_integer(): 

783 return str(int(n)) 

784 return str(n) 

785 

786 

787class SurveyWidePropertyMapPlot(PlotAction): 

788 plotName = pexConfig.Field[str](doc="The name for the plotting task.", optional=True) 

789 

790 def __call__( 

791 self, 

792 data: KeyedData, 

793 plotConfig: SurveyWidePropertyMapAnalysisConfig, 

794 plotInfo: Mapping[str, Union[Mapping[str, str], str, int]], 

795 **kwargs,