LSST Applications g063fba187b+fee0456c91,g0f08755f38+ea96e5a5a3,g1653933729+a8ce1bb630,g168dd56ebc+a8ce1bb630,g1a2382251a+90257ff92a,g20f6ffc8e0+ea96e5a5a3,g217e2c1bcf+937a289c59,g28da252d5a+daa7da44eb,g2bbee38e9b+253935c60e,g2bc492864f+253935c60e,g3156d2b45e+6e55a43351,g32e5bea42b+31359a2a7a,g347aa1857d+253935c60e,g35bb328faa+a8ce1bb630,g3a166c0a6a+253935c60e,g3b1af351f3+a8ce1bb630,g3e281a1b8c+c5dd892a6c,g414038480c+416496e02f,g41af890bb2+afe91b1188,g599934f4f4+0db33f7991,g7af13505b9+e36de7bce6,g80478fca09+da231ba887,g82479be7b0+a4516e59e3,g858d7b2824+ea96e5a5a3,g89c8672015+f4add4ffd5,g9125e01d80+a8ce1bb630,ga5288a1d22+bc6ab8dfbd,gb58c049af0+d64f4d3760,gc28159a63d+253935c60e,gcab2d0539d+3f2b72788c,gcf0d15dbbd+4ea9c45075,gda6a2b7d83+4ea9c45075,gdaeeff99f8+1711a396fd,ge79ae78c31+253935c60e,gef2f8181fd+3031e3cf99,gf0baf85859+c1f95f4921,gfa517265be+ea96e5a5a3,gfa999e8aa5+17cd334064,w.2024.50
LSST Data Management Base Package
Loading...
Searching...
No Matches
Functions | Variables
showVisitSkyMap Namespace Reference

Functions

 bboxToRaDec (bbox, wcs)
 
 getValueAtPercentile (values, percentile=0.5)
 
 get_cmap (n, name="hsv")
 
 main (repo, collections, skymapName=None, tracts=None, visits=None, physicalFilters=None, bands=None, ccds=None, ccdKey="detector", showPatch=False, saveFile=None, showCcds=False, visitVetoFile=None, minOverlapFraction=None, trimToTracts=False, doUnscaledLimitRatio=False, forceScaledLimitRatio=False)
 
 makeWhereInStr (parameterName, parameterList, parameterType)
 
 getTractLimitsDict (skymap, tractList)
 
 getMinMaxLimits (limitsDict)
 
 derivePlotLimits (limitsDict, raToDecLimitRatio=1.0, buffFrac=0.0)
 
 setLimitsToEqualRatio (xMin, xMax, yMin, yMax)
 
 getDetRaDecCorners (ccdKey, ccdId, visit, visitSummary=None, butler=None, doLogWarn=True)
 
 getBand (visitSummary=None, butler=None, visit=None)
 

Variables

 logger = logging.getLogger("lsst.skymap.bin.showVisitSkyMap")
 
 parser = argparse.ArgumentParser()
 
 type
 
 help
 
 str
 
 nargs
 
 metavar
 
 required
 
 default
 
 None
 
 int
 
 action
 
 float
 
 args = parser.parse_args()
 
 level
 
 collections
 
 skymapName
 
 tracts
 
 visits
 
 physicalFilters
 
 bands
 
 ccds
 
 ccdKey
 
 showPatch
 
 saveFile
 
 showCcds
 
 visitVetoFile
 
 minOverlapFraction
 
 trimToTracts
 
 doUnscaledLimitRatio
 
 forceScaledLimitRatio
 

Function Documentation

◆ bboxToRaDec()

showVisitSkyMap.bboxToRaDec ( bbox,
wcs )
Get the corners of a BBox and convert them to lists of RA and Dec.

Definition at line 42 of file showVisitSkyMap.py.

42def bboxToRaDec(bbox, wcs):
43 """Get the corners of a BBox and convert them to lists of RA and Dec."""
44 sphPoints = wcs.pixelToSky(geom.Box2D(bbox).getCorners())
45 ra = [float(sph.getRa().asDegrees()) for sph in sphPoints]
46 dec = [float(sph.getDec().asDegrees()) for sph in sphPoints]
47 return ra, dec
48
49
A floating-point coordinate rectangle geometry.
Definition Box.h:413

◆ derivePlotLimits()

showVisitSkyMap.derivePlotLimits ( limitsDict,
raToDecLimitRatio = 1.0,
buffFrac = 0.0 )
Derive the axis limits to encompass all points in limitsDict.

Parameters
----------
limitsDict : `dict` [`dict`]
    A dictionary keyed on any id. Each entry includes a `dict`
    keyed on "ras" and "decs" including (at least the minimum
    and maximum) RA and Dec values in units of degrees.
raToDecLimitRatio : `float`, optional
    The aspect ratio between RA and Dec to set the plot limits to.  This
    is to namely to set this ratio to that of the focal plane (i.e. such
    that a square detector appears as a square), but any aspect ratio can,
    in principle, be requested.

Returns
-------
xlim, ylim : `tuple` [`float`]
    Two tuples containing the derived min and max values for the x and
    y-axis limits (in degrees), respectively.

Definition at line 614 of file showVisitSkyMap.py.

614def derivePlotLimits(limitsDict, raToDecLimitRatio=1.0, buffFrac=0.0):
615 """Derive the axis limits to encompass all points in limitsDict.
616
617 Parameters
618 ----------
619 limitsDict : `dict` [`dict`]
620 A dictionary keyed on any id. Each entry includes a `dict`
621 keyed on "ras" and "decs" including (at least the minimum
622 and maximum) RA and Dec values in units of degrees.
623 raToDecLimitRatio : `float`, optional
624 The aspect ratio between RA and Dec to set the plot limits to. This
625 is to namely to set this ratio to that of the focal plane (i.e. such
626 that a square detector appears as a square), but any aspect ratio can,
627 in principle, be requested.
628
629 Returns
630 -------
631 xlim, ylim : `tuple` [`float`]
632 Two tuples containing the derived min and max values for the x and
633 y-axis limits (in degrees), respectively.
634 """
635 xLimMin, xLimMax, yLimMin, yLimMax = getMinMaxLimits(limitsDict)
636
637 xDelta0 = xLimMax - xLimMin
638 yDelta0 = yLimMax - yLimMin
639 if raToDecLimitRatio is None:
640 padFrac = 0.05
641 xlim = xLimMax + padFrac*xDelta0, xLimMin - padFrac*xDelta0
642 ylim = yLimMin - padFrac*yDelta0, yLimMax + padFrac*yDelta0
643 return xlim, ylim
644
645 if raToDecLimitRatio == 1.0:
646 if xDelta0 > yDelta0:
647 xLimMin -= buffFrac*yDelta0
648 xLimMax += buffFrac*yDelta0
649 else:
650 yLimMin -= buffFrac*yDelta0
651 yLimMax += buffFrac*yDelta0
652 xLimMin, xLimMax, yLimMin, yLimMax = setLimitsToEqualRatio(xLimMin, xLimMax, yLimMin, yLimMax)
653 else:
654 xLimMin -= buffFrac*xDelta0
655 xLimMax += buffFrac*xDelta0
656 yLimMin -= buffFrac*yDelta0
657 yLimMax += buffFrac*yDelta0
658 xLimMin, xLimMax, yLimMin, yLimMax = setLimitsToEqualRatio(xLimMin, xLimMax, yLimMin, yLimMax)
659 xDelta = xLimMax - xLimMin
660 yDelta = yLimMax - yLimMin
661 if raToDecLimitRatio > 1.0:
662 if yDelta0 > xDelta:
663 xMid = xLimMin + 0.5*(xDelta)
664 xLimMin = xMid - 0.5*yDelta*raToDecLimitRatio
665 xLimMax = xMid + 0.5*yDelta*raToDecLimitRatio
666 else:
667 yMid = yLimMin + 0.5*(yDelta)
668 yLimMin = yMid - 0.5*xDelta/raToDecLimitRatio
669 yLimMax = yMid + 0.5*xDelta/raToDecLimitRatio
670 else:
671 if xDelta0 > yDelta0:
672 yMid = yLimMin + 0.5*(yDelta)
673 yLimMin = yMid - 0.5*xDelta/raToDecLimitRatio
674 yLimMax = yMid + 0.5*xDelta/raToDecLimitRatio
675 else:
676 xMid = xLimMin + 0.5*(xDelta)
677 xLimMin = xMid - 0.5*yDelta*raToDecLimitRatio
678 xLimMax = xMid + 0.5*yDelta*raToDecLimitRatio
679 xlim = xLimMax, xLimMin
680 ylim = yLimMin, yLimMax
681 return xlim, ylim
682
683

◆ get_cmap()

showVisitSkyMap.get_cmap ( n,
name = "hsv" )
Returns a function that maps each index in 0, 1, ..., n-1 to a distinct
RGB color; the keyword argument name must be a standard mpl colormap name.

Definition at line 72 of file showVisitSkyMap.py.

72def get_cmap(n, name="hsv"):
73 """Returns a function that maps each index in 0, 1, ..., n-1 to a distinct
74 RGB color; the keyword argument name must be a standard mpl colormap name.
75 """
76 return matplotlib.colormaps[name].resampled(n)
77
78

◆ getBand()

showVisitSkyMap.getBand ( visitSummary = None,
butler = None,
visit = None )
Determine band and physical filter for given visit.

Parameters
----------
visitSummary : `lsst.afw.table.ExposureCatalog` or `None`, optional
    The visitSummary table for the visit for which to determine the band.
butler : `lsst.daf.butler.Butler` or `None`, optional
    The butler from which to look up the Dimension Records. Only needed
    if ``visitSummary`` is `None`.
visit : `int` or `None, optional
    The visit number for which to determine the band. Only needed
    if ``visitSummary`` is `None`.

Returns
-------
band, physicalFilter : `str`
    The band and physical filter for the given visit.

Definition at line 738 of file showVisitSkyMap.py.

738def getBand(visitSummary=None, butler=None, visit=None):
739 """Determine band and physical filter for given visit.
740
741 Parameters
742 ----------
743 visitSummary : `lsst.afw.table.ExposureCatalog` or `None`, optional
744 The visitSummary table for the visit for which to determine the band.
745 butler : `lsst.daf.butler.Butler` or `None`, optional
746 The butler from which to look up the Dimension Records. Only needed
747 if ``visitSummary`` is `None`.
748 visit : `int` or `None, optional
749 The visit number for which to determine the band. Only needed
750 if ``visitSummary`` is `None`.
751
752 Returns
753 -------
754 band, physicalFilter : `str`
755 The band and physical filter for the given visit.
756 """
757 if visitSummary is not None:
758 band = visitSummary[0]["band"]
759 physicalFilter = visitSummary[0]["physical_filter"]
760 else:
761 record = list(butler.registry.queryDimensionRecords("band", visit=visit))[0]
762 band = record.name
763 record = list(butler.registry.queryDimensionRecords("physical_filter", visit=visit))[0]
764 physicalFilter = record.name
765 return band, physicalFilter
766
767

◆ getDetRaDecCorners()

showVisitSkyMap.getDetRaDecCorners ( ccdKey,
ccdId,
visit,
visitSummary = None,
butler = None,
doLogWarn = True )
Compute the RA/Dec corners lists for a given detector in a visit.

Definition at line 713 of file showVisitSkyMap.py.

713def getDetRaDecCorners(ccdKey, ccdId, visit, visitSummary=None, butler=None, doLogWarn=True):
714 """Compute the RA/Dec corners lists for a given detector in a visit.
715 """
716 raCorners, decCorners = None, None
717 if visitSummary is not None:
718 row = visitSummary.find(ccdId)
719 if row is None:
720 if doLogWarn:
721 logger.warn("No row found for %d in visitSummary of visit %d. "
722 "Skipping and continuing...", ccdId, visit)
723 else:
724 raCorners = list(row["raCorners"])
725 decCorners = list(row["decCorners"])
726 else:
727 try:
728 dataId = {"visit": visit, ccdKey: ccdId}
729 wcs = butler.get("calexp.wcs", dataId)
730 bbox = butler.get("calexp.bbox", dataId)
731 raCorners, decCorners = bboxToRaDec(bbox, wcs)
732 except LookupError as e:
733 logger.warn("%s Skipping and continuing...", e)
734
735 return raCorners, decCorners
736
737

◆ getMinMaxLimits()

showVisitSkyMap.getMinMaxLimits ( limitsDict)
Derive the min and max axis limits of points in limitsDict.

Parameters
----------
limitsDict : `dict` [`dict`]
    A dictionary keyed on any id. Each entry includes a `dict`
    keyed on "ras" and "decs" including (at least the minimum
    and maximum) RA and Dec values in units of degrees.

Returns
-------
xLimMin, xLimMax, yLimMin, yLimMax : `tuple` [`float`]
    The min and max values for the x and y-axis limits, respectively.

Definition at line 588 of file showVisitSkyMap.py.

588def getMinMaxLimits(limitsDict):
589 """Derive the min and max axis limits of points in limitsDict.
590
591 Parameters
592 ----------
593 limitsDict : `dict` [`dict`]
594 A dictionary keyed on any id. Each entry includes a `dict`
595 keyed on "ras" and "decs" including (at least the minimum
596 and maximum) RA and Dec values in units of degrees.
597
598 Returns
599 -------
600 xLimMin, xLimMax, yLimMin, yLimMax : `tuple` [`float`]
601 The min and max values for the x and y-axis limits, respectively.
602 """
603 xLimMin, yLimMin = 1e12, 1e12
604 xLimMax, yLimMax = -1e12, -1e12
605 for limitId, limits in limitsDict.items():
606 xLimMin = min(xLimMin, min(limits["ras"]))
607 xLimMax = max(xLimMax, max(limits["ras"]))
608 yLimMin = min(yLimMin, min(limits["decs"]))
609 yLimMax = max(yLimMax, max(limits["decs"]))
610
611 return xLimMin, xLimMax, yLimMin, yLimMax
612
613
int min
int max

◆ getTractLimitsDict()

showVisitSkyMap.getTractLimitsDict ( skymap,
tractList )
Return a dict containing tract limits needed for outline plotting.

Parameters
----------
skymap : `lsst.skymap.BaseSkyMap`
    The sky map used for this dataset. Used to obtain tract
    parameters.
tractList : `list` [`int`]
    The list of tract ids (as integers) for which to determine the
    limits.

Returns
-------
tractLimitsDict : `dict` [`dict`]
    A dictionary keyed on tract id. Each entry includes a `dict`
    including the tract RA corners, Dec corners, and the tract center,
    all in units of degrees. These are used for plotting the tract
    outlines.

Definition at line 550 of file showVisitSkyMap.py.

550def getTractLimitsDict(skymap, tractList):
551 """Return a dict containing tract limits needed for outline plotting.
552
553 Parameters
554 ----------
555 skymap : `lsst.skymap.BaseSkyMap`
556 The sky map used for this dataset. Used to obtain tract
557 parameters.
558 tractList : `list` [`int`]
559 The list of tract ids (as integers) for which to determine the
560 limits.
561
562 Returns
563 -------
564 tractLimitsDict : `dict` [`dict`]
565 A dictionary keyed on tract id. Each entry includes a `dict`
566 including the tract RA corners, Dec corners, and the tract center,
567 all in units of degrees. These are used for plotting the tract
568 outlines.
569 """
570 tractLimitsDict = {}
571 for tract in tractList:
572 tractInfo = skymap[tract]
573 tractBbox = tractInfo.outer_sky_polygon.getBoundingBox()
574 tractCenter = tractBbox.getCenter()
575 tractRa0 = (tractCenter[0] - tractBbox.getWidth() / 2).asDegrees()
576 tractRa1 = (tractCenter[0] + tractBbox.getWidth() / 2).asDegrees()
577 tractDec0 = (tractCenter[1] - tractBbox.getHeight() / 2).asDegrees()
578 tractDec1 = (tractCenter[1] + tractBbox.getHeight() / 2).asDegrees()
579 tractLimitsDict[tract] = {
580 "ras": [tractRa0, tractRa1, tractRa1, tractRa0, tractRa0],
581 "decs": [tractDec0, tractDec0, tractDec1, tractDec1, tractDec0],
582 "center": [tractCenter[0].asDegrees(), tractCenter[1].asDegrees()],
583 }
584
585 return tractLimitsDict
586
587

◆ getValueAtPercentile()

showVisitSkyMap.getValueAtPercentile ( values,
percentile = 0.5 )
Return a value a fraction of the way between the min and max values in a
list.

Parameters
----------
values : `list` [`float`]
    The list of values under consideration.
percentile : `float`, optional
    The percentile (expressed as a number between 0.0 and 1.0) at which
    to determine the value in `values`.

Returns
-------
result : `float`
    The value at the given percentile of ``values``.

Definition at line 50 of file showVisitSkyMap.py.

50def getValueAtPercentile(values, percentile=0.5):
51 """Return a value a fraction of the way between the min and max values in a
52 list.
53
54 Parameters
55 ----------
56 values : `list` [`float`]
57 The list of values under consideration.
58 percentile : `float`, optional
59 The percentile (expressed as a number between 0.0 and 1.0) at which
60 to determine the value in `values`.
61
62 Returns
63 -------
64 result : `float`
65 The value at the given percentile of ``values``.
66 """
67 m = min(values)
68 interval = max(values) - m
69 return m + percentile*interval
70
71

◆ main()

showVisitSkyMap.main ( repo,
collections,
skymapName = None,
tracts = None,
visits = None,
physicalFilters = None,
bands = None,
ccds = None,
ccdKey = "detector",
showPatch = False,
saveFile = None,
showCcds = False,
visitVetoFile = None,
minOverlapFraction = None,
trimToTracts = False,
doUnscaledLimitRatio = False,
forceScaledLimitRatio = False )

Definition at line 79 of file showVisitSkyMap.py.

82 forceScaledLimitRatio=False):
83 if minOverlapFraction is not None and tracts is None:
84 raise RuntimeError("Must specify --tracts if --minOverlapFraction is set")
85 logger.info("Making butler for collections = %s in repo %s", collections, repo)
86 butler = dafButler.Butler(repo, collections=collections)
87 instrument = butler.find_dataset("camera").dataId["instrument"]
88 detectorSkipList = []
89 # Make a guess at the skymapName if not provided
90 if skymapName is None:
91 if instrument == "HSC":
92 skymapName = "hsc_rings_v1"
93 detectorSkipList = [9] # detector 9 has long been dead for HSC
94 elif instrument == "LSSTCam-imSim":
95 skymapName = "DC2"
96 elif instrument == "LSSTComCamSim":
97 skymapName = "ops_rehersal_prep_2k_v1"
98 elif instrument == "LATISS":
99 skymapName = "latiss_v1"
100 elif instrument == "DECam":
101 skymapName = "decam_rings_v1"
102 else:
103 logger.error("Unknown skymapName for instrument: %s. Must specify --skymapName on command line.",
104 instrument)
105 logger.info("instrument = %s skymapName = %s", instrument, skymapName)
106 camera = butler.get("camera", instrument=instrument)
107 skymap = butler.get("skyMap", instrument=instrument, skymap=skymapName)
108
109 whereStr = ""
110 if tracts is not None:
111 tractStr = makeWhereInStr("tract", tracts, int)
112 whereStr += tractStr
113
114 if physicalFilters is not None:
115 physicalFilterStr = makeWhereInStr("physical_filter", physicalFilters, str)
116 whereStr += " AND " + physicalFilterStr if len(whereStr) else " " + physicalFilterStr
117
118 if bands is not None:
119 bandStr = makeWhereInStr("band", bands, str)
120 whereStr += " AND " + bandStr if len(whereStr) else " " + bandStr
121
122 if len(whereStr) > 1:
123 whereStr = "instrument=\'" + instrument + "\' AND skymap=\'" + skymapName + "\' AND " + whereStr
124 else:
125 whereStr = "instrument=\'" + instrument + "\' AND skymap=\'" + skymapName + "\'"
126 logger.info("Applying the following where clause in dataId search: %s", whereStr)
127
128 visitVetoList = []
129 if visitVetoFile is not None:
130 with open(visitVetoFile) as f:
131 content = f.readlines()
132 visitVetoList = [int(visit.strip()) for visit in content]
133
134 if visits is None:
135 dataRefs = list(butler.registry.queryDatasets("calexp", where=whereStr).expanded())
136 visits = []
137 for dataRef in dataRefs:
138 visit = dataRef.dataId.visit.id
139 if visit not in visits and visit not in visitVetoList:
140 visits.append(visit)
141 visits.sort()
142 logger.info(
143 "List of visits (N=%d) satisfying where and veto clauses: %s", len(visits), visits
144 )
145 else:
146 if len(visitVetoList) > 1:
147 visitListTemp = visits.copy()
148 for visit in visitListTemp:
149 if visit in visitVetoList:
150 visits.remove(visit)
151 logger.info("List of visits (N=%d) excluding veto list: %s}", len(visits), visits)
152 logger.info("List of visits (N=%d): %s", len(visits), visits)
153
154 ccdIdList = []
155 for ccd in camera:
156 ccdId = ccd.getId()
157 if ((ccds is None or ccdId in ccds) and ccd.getType() == cameraGeom.DetectorType.SCIENCE
158 and ccdId not in detectorSkipList):
159 ccdIdList.append(ccdId)
160 ccdIdList.sort()
161 nDetTot = len(ccdIdList)
162
163 visitIncludeList = []
164 # Determine the fraction of detectors that overlap any tract under
165 # consideration. If this fraction does not exceed minOverlapFraction,
166 # skip this visit.
167 if minOverlapFraction is not None:
168 for i_v, visit in enumerate(visits):
169 ccdOverlapList = []
170 try:
171 visitSummary = butler.get("visitSummary", visit=visit)
172 except LookupError as e:
173 logger.warn("%s Will try to get wcs from calexp.", e)
174 visitSummary = None
175 if tracts is not None:
176 for tract in tracts:
177 tractInfo = skymap[tract]
178 sphCorners = tractInfo.wcs.pixelToSky(geom.Box2D(tractInfo.bbox).getCorners())
179 tractConvexHull = sphgeom.ConvexPolygon.convexHull(
180 [coord.getVector() for coord in sphCorners])
181 for ccdId in ccdIdList:
182 if ccdId not in ccdOverlapList:
183 raCorners, decCorners = getDetRaDecCorners(
184 ccdKey, ccdId, visit, visitSummary=visitSummary, butler=butler,
185 doLogWarn=False)
186 if raCorners is not None and decCorners is not None:
187 detSphCorners = []
188 for ra, dec in zip(raCorners, decCorners):
189 pt = geom.SpherePoint(geom.Angle(ra, geom.degrees),
190 geom.Angle(dec, geom.degrees))
191 detSphCorners.append(pt)
192 detConvexHull = sphgeom.ConvexPolygon.convexHull(
193 [coord.getVector() for coord in detSphCorners])
194 if tractConvexHull.contains(detConvexHull):
195 ccdOverlapList.append(ccdId)
196
197 if len(ccdOverlapList)/nDetTot >= minOverlapFraction:
198 break
199 if len(ccdOverlapList)/nDetTot < minOverlapFraction:
200 logger.info("Fraction of detectors overlaping any tract for visit %d (%.2f) < "
201 "minimum required (%.2f). Skipping visit...",
202 visit, len(ccdOverlapList)/nDetTot, minOverlapFraction)
203 else:
204 if visit not in visitIncludeList:
205 visitIncludeList.append(visit)
206 else:
207 visitIncludeList = visits
208
209 # Draw the CCDs.
210 ras, decs = [], []
211 bboxesPlotted = []
212 cmap = get_cmap(len(visitIncludeList))
213 alphaEdge = 0.7
214 maxVisitForLegend = 20
215 finalVisitList = []
216 includedBands = []
217 includedPhysicalFilters = []
218 for i_v, visit in enumerate(visitIncludeList):
219 print("Working on visit %d [%d of %d]" % (visit, i_v + 1, len(visitIncludeList)), end="\r")
220 inLegend = False
221 color = cmap(i_v)
222 fillKwargs = {"fill": False, "alpha": alphaEdge, "facecolor": None, "edgecolor": color, "lw": 0.6}
223 try:
224 visitSummary = butler.get("visitSummary", visit=visit)
225 except Exception as e:
226 logger.warn("%s Will try to get wcs from calexp.", e)
227 visitSummary = None
228
229 band, physicalFilter = getBand(visitSummary=visitSummary, butler=butler, visit=visit)
230 if band not in includedBands:
231 includedBands.append(band)
232 if physicalFilter not in includedPhysicalFilters:
233 includedPhysicalFilters.append(physicalFilter)
234
235 for ccdId in ccdIdList:
236 raCorners, decCorners = getDetRaDecCorners(
237 ccdKey, ccdId, visit, visitSummary=visitSummary, butler=butler)
238 if raCorners is not None and decCorners is not None:
239 ras += raCorners
240 decs += decCorners
241 if not inLegend and len(visitIncludeList) <= maxVisitForLegend:
242 plt.fill(raCorners, decCorners, label=str(visit), **fillKwargs)
243 inLegend = True
244 else:
245 plt.fill(raCorners, decCorners, **fillKwargs)
246 plt.fill(raCorners, decCorners, fill=True, alpha=alphaEdge/4, color=color,
247 edgecolor=color)
248 if visit not in finalVisitList:
249 finalVisitList.append(visit)
250 # add CCD serial numbers
251 if showCcds:
252 overlapFrac = 0.2
253 deltaRa = max(raCorners) - min(raCorners)
254 deltaDec = max(decCorners) - min(decCorners)
255 minPoint = geom.Point2D(min(raCorners) + overlapFrac*deltaRa,
256 min(decCorners) + overlapFrac*deltaDec)
257 maxPoint = geom.Point2D(max(raCorners) - overlapFrac*deltaRa,
258 max(decCorners) - overlapFrac*deltaDec)
259 # Use doubles in Box2D to check overlap
260 bboxDouble = geom.Box2D(minPoint, maxPoint)
261 overlaps = [not bboxDouble.overlaps(otherBbox) for otherBbox in bboxesPlotted]
262 if all(overlaps):
263 plt.text(getValueAtPercentile(raCorners), getValueAtPercentile(decCorners),
264 str(ccdId), fontsize=6, ha="center", va="center", color="darkblue")
265 bboxesPlotted.append(bboxDouble)
266
267 logger.info("Final list of visits (N=%d) satisfying where and minOverlapFraction clauses: %s",
268 len(finalVisitList), finalVisitList)
269
270 raToDecLimitRatio = None
271 if len(ras) > 0:
272 tractList = list(set(skymap.findTractIdArray(ras, decs, degrees=True)))
273 minVisitRa, maxVisitRa = min(ras), max(ras)
274 minVisitDec, maxVisitDec = min(decs), max(decs)
275 raVisitDiff = maxVisitRa - minVisitRa
276 decVisitDiff = maxVisitDec - minVisitDec
277 midVisitRa = minVisitRa + 0.5*raVisitDiff
278 midVisitDec = minVisitDec + 0.5*decVisitDiff
279 midRa = np.atleast_1d((midVisitRa*units.deg).to(units.radian).value).astype(np.float64)
280 midDec = np.atleast_1d((midVisitDec*units.deg).to(units.radian).value).astype(np.float64)
281 midSkyCoord = SkyCoord(midVisitRa*units.deg, midVisitDec*units.deg)
282 else:
283 if tracts is not None:
284 logger.info("No calexps were found, but --tracts list was provided, so will go ahead and "
285 "plot the empty tracts.")
286 tractList = tracts
287 trimToTracts = True
288 else:
289 raise RuntimeError("No data to plot (if you want to plot empty tracts, include them as "
290 "a blank-space separated list to the --tracts option.")
291 tractList.sort()
292 logger.info("List of tracts overlapping data: %s", tractList)
293 tractLimitsDict = getTractLimitsDict(skymap, tractList)
294
295 if forceScaledLimitRatio:
296 doUnscaledLimitRatio = False
297 else:
298 # Roughly compute radius in degrees of a single detector. If RA/Dec
299 # coverage is more than 30 times the detector radius, and the RA/Dec
300 # limit ratio is greater than raDecScaleThresh, don't try to scale to
301 # detector coords.
302 radiusMm = camera.computeMaxFocalPlaneRadius()
303 fpRadiusPt = geom.Point2D(radiusMm, radiusMm)
304 focalPlaneToFieldAngle = camera.getTransformMap().getTransform(
305 cameraGeom.FOCAL_PLANE, cameraGeom.FIELD_ANGLE
306 )
307 fpRadiusDeg = np.rad2deg(focalPlaneToFieldAngle.applyForward(fpRadiusPt))[0]
308 detectorRadiusDeg = fpRadiusDeg/np.sqrt(len(camera))
309
310 if trimToTracts:
311 xLimMin, xLimMax, yLimMin, yLimMax = getMinMaxLimits(tractLimitsDict)
312 xDelta0 = xLimMax - xLimMin
313 yDelta0 = yLimMax - yLimMin
314 else:
315 xDelta0 = raVisitDiff
316 yDelta0 = decVisitDiff
317 yLimMin = minVisitDec
318 yLimMax = maxVisitDec
319 raDecScaleThresh = 1.5 # This is a best guess with current testing.
320 if (
321 (xDelta0/yDelta0 > raDecScaleThresh or yDelta0/xDelta0 > raDecScaleThresh)
322 and max(xDelta0, yDelta0) > 70*detectorRadiusDeg
323 and yLimMin < 75.0 and yLimMax > -75.0
324 ):
325 logger.info(
326 "Sky coverage is large (and not too close to a pole), so not scaling to detector coords."
327 )
328 doUnscaledLimitRatio = True
329
330 if not doUnscaledLimitRatio:
331 # Find a detector that contains the mid point in RA/Dec (or the closest
332 # one) to set the plot aspect ratio.
333 minDistToMidCoord = 1e12
334 minSepVisit = None
335 minSepCcdId = None
336 for i_v, visit in enumerate(visits):
337 try:
338 visitSummary = butler.get("visitSummary", visit=visit)
339 except Exception as e:
340 logger.warn("%s Will try to get wcs from calexp.", e)
341 visitSummary = None
342 for ccdId in ccdIdList:
343 raCorners, decCorners = getDetRaDecCorners(
344 ccdKey, ccdId, visit, visitSummary=visitSummary, butler=butler, doLogWarn=False)
345 if raCorners is not None and decCorners is not None:
346 detSphCorners = []
347 for ra, dec in zip(raCorners, decCorners):
348 pt = geom.SpherePoint(geom.Angle(ra, geom.degrees),
349 geom.Angle(dec, geom.degrees))
350 detSphCorners.append(pt)
351 ptSkyCoord = SkyCoord(ra*units.deg, dec*units.deg)
352 separation = (midSkyCoord.separation(ptSkyCoord)).degree
353 if separation < minDistToMidCoord:
354 minSepVisit = visit
355 minSepCcdId = ccdId
356 minDistToMidCoord = separation
357 detConvexHull = sphgeom.ConvexPolygon(
358 [coord.getVector() for coord in detSphCorners])
359 if detConvexHull.contains(midRa, midDec) and raToDecLimitRatio is None:
360 logger.info(
361 "visit/det overlapping plot coord mid point in RA/Dec: %d %d", visit, ccdId
362 )
363 raToDecLimitRatio = (
364 (max(raCorners) - min(raCorners))/(max(decCorners) - min(decCorners))
365 )
366 det = camera[ccdId]
367 width = det.getBBox().getWidth()
368 height = det.getBBox().getHeight()
369 if raToDecLimitRatio > 1.0:
370 raToDecLimitRatio /= max(height/width, width/height)
371 else:
372 if raToDecLimitRatio < 1.0:
373 raToDecLimitRatio *= max(height/width, width/height)
374 break
375 if raToDecLimitRatio is not None:
376 break
377
378 if raToDecLimitRatio is None and minSepVisit is not None:
379 try:
380 visitSummary = butler.get("visitSummary", visit=minSepVisit)
381 except Exception as e:
382 logger.warn("%s Will try to get wcs from calexp.", e)
383 visitSummary = None
384 raCorners, decCorners = getDetRaDecCorners(
385 ccdKey, minSepCcdId, minSepVisit, visitSummary=visitSummary, butler=butler, doLogWarn=False)
386 for ra, dec in zip(raCorners, decCorners):
387 pt = geom.SpherePoint(geom.Angle(ra, geom.degrees),
388 geom.Angle(dec, geom.degrees))
389 detSphCorners.append(pt)
390 detConvexHull = sphgeom.ConvexPolygon([coord.getVector() for coord in detSphCorners])
391 logger.info(
392 "visit/det closest to plot coord mid point in RA/Dec (none actually overlap it): %d %d",
393 minSepVisit, minSepCcdId
394 )
395 raToDecLimitRatio = (max(raCorners) - min(raCorners))/(max(decCorners) - min(decCorners))
396 det = camera[minSepCcdId]
397 width = det.getBBox().getWidth()
398 height = det.getBBox().getHeight()
399 if raToDecLimitRatio > 1.0:
400 raToDecLimitRatio /= max(height/width, width/height)
401 else:
402 if raToDecLimitRatio < 1.0:
403 raToDecLimitRatio *= max(height/width, width/height)
404
405 if trimToTracts is True:
406 xlim, ylim = derivePlotLimits(tractLimitsDict, raToDecLimitRatio=raToDecLimitRatio, buffFrac=0.04)
407 else:
408 visitLimitsDict = {"allVisits": {"ras": [minVisitRa, maxVisitRa], "decs": [minVisitDec, maxVisitDec]}}
409 xlim, ylim = derivePlotLimits(visitLimitsDict, raToDecLimitRatio=raToDecLimitRatio, buffFrac=0.04)
410
411 if doUnscaledLimitRatio:
412 boxAspectRatio = abs((ylim[1] - ylim[0])/(xlim[1] - xlim[0]))
413 else:
414 boxAspectRatio = 1.0
415
416 # Draw the skymap.
417 alpha0 = 1.0
418 tractHandleList = []
419 tractStrList = []
420 if tracts is not None:
421 tractOutlineList = list(set(tracts + tractList))
422 else:
423 tractOutlineList = tractList
424 tractOutlineList.sort()
425 logger.info("List of tract outlines being plotted: %s", tractOutlineList)
426 for i_t, tract in enumerate(tractOutlineList):
427 alpha = max(0.1, alpha0 - i_t*1.0/len(tractOutlineList))
428 tractInfo = skymap[tract]
429 tCenter = tractInfo.ctr_coord
430 tCenterRa = tCenter.getRa().asDegrees()
431 tCenterDec = tCenter.getDec().asDegrees()
432 fracDeltaX = 0.02*abs((xlim[1] - xlim[0]))
433 fracDeltaY = 0.02*abs((ylim[1] - ylim[0]))
434 if (xlim[1] + fracDeltaX < tCenterRa < xlim[0] - fracDeltaX
435 and ylim[0] + fracDeltaY < tCenterDec < ylim[1] - fracDeltaY):
436 if len(tractOutlineList) > 1 or not showPatch:
437 if not showPatch:
438 plt.text(tCenterRa, tCenterDec, tract, fontsize=7, alpha=alpha, ha="center", va="center")
439 else:
440 plt.text(tCenterRa, tCenterDec, tract, fontsize=7, alpha=1, color="white",
441 path_effects=[pathEffects.withStroke(linewidth=3, foreground="black")],
442 fontweight=500, ha="center", va="center", zorder=5)
443 ra, dec = bboxToRaDec(tractInfo.bbox, tractInfo.getWcs())
444 plt.fill(ra, dec, fill=False, edgecolor="k", lw=1, linestyle="dashed", alpha=alpha)
445 tractArtist = matplotlib.patches.Patch(fill=False, edgecolor="k", linestyle="dashed", alpha=alpha)
446 tractHandleList.append(tractArtist)
447 tractStrList.append(str(tract))
448 if showPatch:
449 patchColor = "k"
450 for patch in tractInfo:
451 ra, dec = bboxToRaDec(patch.getInnerBBox(), tractInfo.getWcs())
452 plt.fill(ra, dec, fill=False, edgecolor=patchColor, lw=0.5, linestyle=(0, (5, 6)),
453 alpha=alpha)
454 if (xlim[1] + fracDeltaX < getValueAtPercentile(ra) < xlim[0] - fracDeltaX
455 and ylim[0] + fracDeltaY < getValueAtPercentile(dec) < ylim[1] - fracDeltaY):
456 plt.text(getValueAtPercentile(ra), getValueAtPercentile(dec),
457 str(patch.sequential_index), fontsize=5, color=patchColor,
458 ha="center", va="center", alpha=alpha)
459
460 # Add labels and save.
461 ax = plt.gca()
462 ax.set_xlim(xlim)
463 ax.set_ylim(ylim)
464 ax.set_box_aspect(boxAspectRatio)
465 if abs(xlim[1] > 99.99):
466 ax.tick_params("x", labelrotation=45, pad=0, labelsize=8)
467 else:
468 ax.tick_params("x", labelrotation=0, pad=0, labelsize=8)
469 ax.tick_params("y", labelsize=8)
470 ax.set_xlabel("RA (deg)", fontsize=9)
471 ax.set_ylabel("Dec (deg)", fontsize=9)
472
473 visitScaleOffset = None
474 if len(visitIncludeList) > maxVisitForLegend:
475 nz = matplotlib.colors.Normalize()
476 colorBarScale = finalVisitList
477 if max(finalVisitList) > 9999999:
478 visitScaleOffset = min(finalVisitList)
479 colorBarScale = [visit - visitScaleOffset for visit in finalVisitList]
480 nz.autoscale(colorBarScale)
481 cax, _ = matplotlib.colorbar.make_axes(plt.gca(), pad=0.03)
482 cax.tick_params(labelsize=7)
483 cb = matplotlib.colorbar.ColorbarBase(cax, cmap=cmap, norm=nz, alpha=alphaEdge,
484 format=lambda x, _: f"{x:.0f}")
485 cb.ax.yaxis.get_offset_text().set_fontsize(7)
486 colorBarLabel = "visit number"
487 if visitScaleOffset is not None:
488 colorBarLabel += " - {:d}".format(visitScaleOffset)
489 cb.set_label(colorBarLabel, rotation=-90, labelpad=13, fontsize=9)
490 tractLegend = Legend(ax, tractHandleList, tractStrList, loc="upper right", fancybox=True,
491 shadow=True, fontsize=5, title_fontsize=6, title="tracts")
492 ax.add_artist(tractLegend)
493 else:
494 if len(visitIncludeList) > 0:
495 xBboxAnchor = min(1.25, max(1.03, boxAspectRatio*1.15))
496 ax.legend(loc="center left", bbox_to_anchor=(xBboxAnchor, 0.5), fancybox=True,
497 shadow=True, fontsize=6, title_fontsize=6, title="visits")
498 # Create the second legend and add the artist manually.
499 tractLegend = Legend(ax, tractHandleList, tractStrList, loc="center left", bbox_to_anchor=(1.0, 0.5),
500 fancybox=True, shadow=True, fontsize=6, title_fontsize=6, title="tracts")
501 ax.add_artist(tractLegend)
502
503 titleStr = repo + "\n" + collections[0]
504 if len(collections) > 1:
505 for collection in collections[1:]:
506 titleStr += "\n" + collection
507 titleStr += "\nnVisit: {}".format(str(len(finalVisitList)))
508 if minOverlapFraction is not None:
509 titleStr += " (minOvlpFrac = {:.2f})".format(minOverlapFraction)
510 if len(includedBands) > 0:
511 titleStr += " bands: {}".format(str(includedBands).translate({ord(i): None for i in "[]\'"}))
512 if len(includedPhysicalFilters) > 0:
513 if len(includedPhysicalFilters[0]) > 9:
514 titleStr += "\n"
515 titleStr += " physical filters: {}".format(str(includedPhysicalFilters).translate(
516 {ord(i): None for i in "[]\'"}))
517 ax.set_title("{}".format(titleStr), fontsize=8)
518
519 fig = plt.gcf()
520 if boxAspectRatio > 1.0:
521 minInches = max(4.0, 0.3*abs(xlim[1] - xlim[0]))
522 xInches = minInches
523 yInches = min(120.0, boxAspectRatio*minInches)
524 fig.set_size_inches(xInches, yInches)
525 if boxAspectRatio < 1.0:
526 minInches = max(4.0, 0.3*abs(ylim[1] - ylim[0]))
527 xInches = min(120.0, minInches/boxAspectRatio)
528 yInches = minInches
529 fig.set_size_inches(xInches, yInches)
530 if saveFile is not None:
531 logger.info("Saving file in: %s", saveFile)
532 fig.savefig(saveFile, bbox_inches="tight", dpi=150)
533 else:
534 fig.show()
535
536
A class representing an angle.
Definition Angle.h:128
Point in an unspecified spherical coordinate system.
Definition SpherePoint.h:57

◆ makeWhereInStr()

showVisitSkyMap.makeWhereInStr ( parameterName,
parameterList,
parameterType )
Create the string to be used in the where clause for registry lookup.

Definition at line 537 of file showVisitSkyMap.py.

537def makeWhereInStr(parameterName, parameterList, parameterType):
538 """Create the string to be used in the where clause for registry lookup.
539 """
540 typeStr = "\'" if parameterType is str else ""
541 whereInStr = parameterName + " IN (" + typeStr + str(parameterList[0]) + typeStr
542 if len(parameterList) > 1:
543 for param in parameterList[1:]:
544 whereInStr += ", " + typeStr + str(param) + typeStr
545 whereInStr += ")"
546
547 return whereInStr
548
549

◆ setLimitsToEqualRatio()

showVisitSkyMap.setLimitsToEqualRatio ( xMin,
xMax,
yMin,
yMax )
For a given set of x/y min/max, redefine to have equal aspect ratio.

The limits are extended on both ends such that the central value is
preserved.

Parameters
----------
xMin, xMax, yMin, yMax : `float`
    The min/max values of the x/y ranges for which to match in dynamic
    range while perserving the central values.

Returns
-------
xMin, xMax, yMin, yMax : `float`
    The adjusted min/max values of the x/y ranges with equal aspect ratios.

Definition at line 684 of file showVisitSkyMap.py.

684def setLimitsToEqualRatio(xMin, xMax, yMin, yMax):
685 """For a given set of x/y min/max, redefine to have equal aspect ratio.
686
687 The limits are extended on both ends such that the central value is
688 preserved.
689
690 Parameters
691 ----------
692 xMin, xMax, yMin, yMax : `float`
693 The min/max values of the x/y ranges for which to match in dynamic
694 range while perserving the central values.
695
696 Returns
697 -------
698 xMin, xMax, yMin, yMax : `float`
699 The adjusted min/max values of the x/y ranges with equal aspect ratios.
700 """
701 xDelta = xMax - xMin
702 yDelta = yMax - yMin
703 deltaDiff = yDelta - xDelta
704 if deltaDiff > 0:
705 xMin -= 0.5 * deltaDiff
706 xMax += 0.5 * deltaDiff
707 elif deltaDiff < 0:
708 yMin -= 0.5 * abs(deltaDiff)
709 yMax += 0.5 * abs(deltaDiff)
710 return xMin, xMax, yMin, yMax
711
712

Variable Documentation

◆ action

showVisitSkyMap.action

Definition at line 790 of file showVisitSkyMap.py.

◆ args

showVisitSkyMap.args = parser.parse_args()

Definition at line 809 of file showVisitSkyMap.py.

◆ bands

showVisitSkyMap.bands

Definition at line 812 of file showVisitSkyMap.py.

◆ ccdKey

showVisitSkyMap.ccdKey

Definition at line 812 of file showVisitSkyMap.py.

◆ ccds

showVisitSkyMap.ccds

Definition at line 812 of file showVisitSkyMap.py.

◆ collections

showVisitSkyMap.collections

Definition at line 811 of file showVisitSkyMap.py.

◆ default

showVisitSkyMap.default

Definition at line 775 of file showVisitSkyMap.py.

◆ doUnscaledLimitRatio

showVisitSkyMap.doUnscaledLimitRatio

Definition at line 815 of file showVisitSkyMap.py.

◆ float

showVisitSkyMap.float

Definition at line 799 of file showVisitSkyMap.py.

◆ forceScaledLimitRatio

showVisitSkyMap.forceScaledLimitRatio

Definition at line 816 of file showVisitSkyMap.py.

◆ help

showVisitSkyMap.help

Definition at line 771 of file showVisitSkyMap.py.

◆ int

showVisitSkyMap.int

Definition at line 776 of file showVisitSkyMap.py.

◆ level

showVisitSkyMap.level

Definition at line 810 of file showVisitSkyMap.py.

◆ logger

showVisitSkyMap.logger = logging.getLogger("lsst.skymap.bin.showVisitSkyMap")

Definition at line 39 of file showVisitSkyMap.py.

◆ metavar

showVisitSkyMap.metavar

Definition at line 774 of file showVisitSkyMap.py.

◆ minOverlapFraction

showVisitSkyMap.minOverlapFraction

Definition at line 814 of file showVisitSkyMap.py.

◆ nargs

showVisitSkyMap.nargs

Definition at line 772 of file showVisitSkyMap.py.

◆ None

showVisitSkyMap.None

Definition at line 775 of file showVisitSkyMap.py.

◆ parser

showVisitSkyMap.parser = argparse.ArgumentParser()

Definition at line 769 of file showVisitSkyMap.py.

◆ physicalFilters

showVisitSkyMap.physicalFilters

Definition at line 812 of file showVisitSkyMap.py.

◆ required

showVisitSkyMap.required

Definition at line 774 of file showVisitSkyMap.py.

◆ saveFile

showVisitSkyMap.saveFile

Definition at line 813 of file showVisitSkyMap.py.

◆ showCcds

showVisitSkyMap.showCcds

Definition at line 813 of file showVisitSkyMap.py.

◆ showPatch

showVisitSkyMap.showPatch

Definition at line 813 of file showVisitSkyMap.py.

◆ skymapName

showVisitSkyMap.skymapName

Definition at line 811 of file showVisitSkyMap.py.

◆ str

showVisitSkyMap.str

Definition at line 772 of file showVisitSkyMap.py.

◆ tracts

showVisitSkyMap.tracts

Definition at line 811 of file showVisitSkyMap.py.

◆ trimToTracts

showVisitSkyMap.trimToTracts

Definition at line 815 of file showVisitSkyMap.py.

◆ type

showVisitSkyMap.type

Definition at line 770 of file showVisitSkyMap.py.

◆ visits

showVisitSkyMap.visits

Definition at line 811 of file showVisitSkyMap.py.

◆ visitVetoFile

showVisitSkyMap.visitVetoFile

Definition at line 814 of file showVisitSkyMap.py.