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