Coverage for python / lsst / analysis / tools / scripts / allSkyProj.py: 0%
70 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-17 09:36 +0000
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/>.
22from typing import Any
24import matplotlib
25import matplotlib.cm as mplcm
26import matplotlib.pyplot as plt
27import numpy as np
28import skyproj # Eli's matplotlib + astronomical projections library
30from lsst.daf.butler import Butler
31from lsst.geom import SpherePoint, degrees
32from lsst.sphgeom import Box, ConvexPolygon, LonLat
35class RegionDisplay:
36 def __init__(self, proj: skyproj.Skyproj):
37 self.proj = proj
38 self.bbox = Box()
40 def add_polygon(self, region: ConvexPolygon, update_bbox: bool = True, **kwargs: Any) -> None:
41 """Plot a sphgeom ConvexPolygon on a sky plot."""
42 vertices = np.array(
43 [
44 (LonLat.longitudeOf(v).asDegrees(), LonLat.latitudeOf(v).asDegrees())
45 for v in region.getVertices()
46 ],
47 dtype=float,
48 )
49 self.proj.fill(vertices[:, 0], vertices[:, 1], **kwargs)
50 if update_bbox:
51 self.bbox.expandTo(region.getBoundingBox())
53 def add_polygon_corners(self, corners: list, **kwargs: Any) -> None:
54 """Plot a sphgeom ConvexPolygon on a sky plot."""
55 corners = np.array(corners)
56 self.proj.fill(corners[:, 0], corners[:, 1], **kwargs)
57 corner_coords = []
58 for corner in corners:
59 corner_coords.append(SpherePoint(corner[0], corner[1], units=degrees).getVector())
60 polygon = ConvexPolygon(corner_coords)
61 self.bbox.expandTo(polygon.getBoundingBox())
63 def update_bounds(self) -> None:
64 """Update the bounds of the plot to just fit all regions."""
65 # print(self.bbox)
66 lon = self.bbox.getLon()
67 lat = self.bbox.getLat()
68 self.proj.set_extent(
69 [lon.getA().asDegrees(), lon.getB().asDegrees(), lat.getA().asDegrees(), lat.getB().asDegrees()]
70 )
73# collection: str = "HSC/runs/RC2/w_2023_07/DM-38042"
74# butlerPath: str = "/repo/main"
75# collection: str = "2.2i/runs/test-med-1/w_2023_13/DM-38781"
76# butlerPath: str = "/repo/dc2"
78collection: str = "u/csaunder/DM-39933/analysis_tools_baseline"
79butlerPath: str = "/sdf/group/rubin/repo/main_20210215"
81butler = Butler(butlerPath)
82print("butler loaded")
84dsTypes = butler.registry.queryDatasetTypes()
85print("dsTypes loaded")
86# print([dsType for dsType in dsTypes])
87accumulator = {}
88# dataID = {'instrument':'HSC', 'skymap':'hsc_rings_v1', 'band':'i',
89# 'visit':1228}
90# m = butler.get('sourceTable_visit_gaia_dr2_20200414_match_metrics',
91# collections=[collection], dataId=dataID)
92# print(m)
93refs = butler.registry.queryDatasets(
94 "sourceTable_visit_gaia_dr2_20200414_match_metrics", collections=[collection], findFirst=True
95)
96for ref in refs:
97 if ref.dataId["band"] == "i":
98 bundle = butler.get(ref)
99 # print(bundle)
100 # print(ref.dataId['band'])
101 dataId = ref.dataId
102 coords = butler.get(
103 "sourceTable_visit",
104 dataId=dataId,
105 collections=[collection],
106 parameters={"columns": ["coord_ra", "coord_dec"]},
107 )
108 min_ra = np.nanmin(coords["coord_ra"])
109 max_ra = np.nanmax(coords["coord_ra"])
110 min_dec = np.nanmin(coords["coord_dec"])
111 max_dec = np.nanmax(coords["coord_dec"])
112 corners = [(min_ra, min_dec), (max_ra, min_dec), (max_ra, max_dec), (min_ra, max_dec)]
113 for name, measurements in bundle.items():
114 for measurement in measurements:
115 compoundName = name
116 # print(name)
117 container = accumulator.setdefault(compoundName, list())
118 container.append((measurement.quantity.value, corners))
119 # print(accumulator)
120"""for dsType in dsTypes:
121 if dsType.storageClass_name != "MetricValue":#"MetricMeasurementBundle":
122 continue
123 if "tract" not in dsType.dimensions:
124 continue
125 refs = list(butler.registry.queryDatasets(dsType, collections=[collection],
126 findFirst=True).expanded())
127 #print(refs)
128"""
129# for ref in refs:
130# print(ref)
131# bundle = butler.get(ref)
132# tmp.append(bundle.metric_name)
133# print('bundle.name: ', bundle.metric_name)
134# help(bundle)
135# for name, measurements in bundle.items():
136# if 'matchedVisitCore' in str(bundle.metric_name):
137# #print(bundle.metric_name)
138# """for measurement in measurements:
139# compoundName = f"{dsType.name}_{name}_{measurement.metric_name.metric}
140# "
141# container = accumulator.setdefault(compoundName, list())
142# container.append((measurement.quantity.value, ref.dataId.region))"""
143# print(set(tmp))
146for name, values in accumulator.items():
147 minimum = min(v for v, _ in values)
148 maximum = max(v for v, _ in values)
149 norm = matplotlib.colors.Normalize(vmin=minimum, vmax=maximum, clip=True)
150 mapper = mplcm.ScalarMappable(norm=norm)
152 fig, ax = plt.subplots(figsize=(8, 5))
153 sp = skyproj.McBrydeSkyproj(ax=ax)
154 regionDisplay = RegionDisplay(sp)
155 for value, region in values:
156 regionDisplay.add_polygon_corners(region, color=mapper.to_rgba(value))
157 regionDisplay.update_bounds()
158 plt.colorbar(mappable=mapper, ax=ax)
159 plt.title(name)
160 # plt.savefig(f"/sdf/home/n/nate2/testplots2/{name}_skyproj.png")
161 plt.savefig(f"/sdf/home/m/mccann/rubin-user/DM-40203/testplots/{name}_skyproj.png")
162 print(str(name) + "_skyproj.png saved")
164 plt.close()