Coverage for python/lsst/summit/utils/utils.py: 16%
622 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-13 02:47 -0800
1# This file is part of summit_utils.
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
23import datetime
24import logging
25import os
26from collections.abc import Iterable, Sequence
27from dataclasses import dataclass
28from typing import TYPE_CHECKING, Any
30import astropy.units as u
31import matplotlib
32import numpy as np
33import numpy.typing as npt
34import statsmodels.api as sm # type: ignore[import-untyped]
35from astro_metadata_translator import ObservationInfo
36from astropy.coordinates import AltAz, EarthLocation, SkyCoord, get_sun
37from astropy.stats import mad_std
38from astropy.time import Time
39from dateutil.tz import gettz
40from deprecated.sphinx import deprecated
41from matplotlib.patches import Rectangle
42from scipy.ndimage import gaussian_filter
43from sklearn.linear_model import LinearRegression, RANSACRegressor # type: ignore[import-untyped]
45import lsst.afw.detection as afwDetect
46import lsst.afw.detection as afwDetection
47import lsst.afw.image as afwImage
48import lsst.afw.math as afwMath
49import lsst.daf.base as dafBase
50import lsst.daf.butler as dafButler
51import lsst.geom as geom
52import lsst.obs.lsst.translators
53import lsst.pipe.base as pipeBase
54import lsst.utils.packages as packageUtils
55from lsst.afw.coord import Weather
56from lsst.afw.detection import Footprint, FootprintSet
57from lsst.daf.butler.cli.cliLog import CliLog
58from lsst.obs.lsst import Latiss, LsstCam, LsstComCam, LsstComCamSim
59from lsst.obs.lsst.translators.latiss import AUXTEL_LOCATION
60from lsst.obs.lsst.translators.lsst import FILTER_DELIMITER, SIMONYI_LOCATION
62from .astrometry.utils import genericCameraHeaderToWcs
64if TYPE_CHECKING:
65 from lsst.afw.cameraGeom import Camera
67__all__ = [
68 "SIGMATOFWHM",
69 "FWHMTOSIGMA",
70 "EFD_CLIENT_MISSING_MSG",
71 "GOOGLE_CLOUD_MISSING_MSG",
72 "AUXTEL_LOCATION",
73 "countPixels",
74 "quickSmooth",
75 "argMax2d",
76 "dayObsIntToString",
77 "dayObsSeqNumToVisitId",
78 "getImageStats",
79 "detectObjectsInExp",
80 "fluxesFromFootprints",
81 "fluxFromFootprint",
82 "humanNameForCelestialObject",
83 "getFocusFromHeader",
84 "checkStackSetup",
85 "setupLogging",
86 "getCurrentDayObs_datetime",
87 "getCurrentDayObs_int",
88 "getCurrentDayObs_humanStr",
89 "getExpRecordAge",
90 "getSite",
91 "getAltAzFromSkyPosition",
92 "getExpPositionOffset",
93 "starTrackerFileToExposure",
94 "obsInfoToDict",
95 "getFieldNameAndTileNumber",
96 "getAirmassSeeingCorrection",
97 "getFilterSeeingCorrection", # deprecated
98 "getBandpassSeeingCorrection",
99 "getCdf",
100 "getQuantiles",
101 "digitizeData",
102 "getBboxAround",
103 "bboxToMatplotlibRectanle",
104 "computeExposureId",
105 "computeCcdExposureId",
106 "getDetectorIds",
107 "getImageArray",
108 "getSunAngle",
109 "RobustFitter",
110]
113SIGMATOFWHM = 2.0 * np.sqrt(2.0 * np.log(2.0))
114FWHMTOSIGMA = 1 / SIGMATOFWHM
116EFD_CLIENT_MISSING_MSG = (
117 "ImportError: lsst_efd_client not found. Please install with:\n" " pip install lsst-efd-client"
118)
120GOOGLE_CLOUD_MISSING_MSG = (
121 "ImportError: Google cloud storage not found. Please install with:\n"
122 " pip install google-cloud-storage"
123)
126def countPixels(maskedImage: afwImage.MaskedImage, maskPlane: str) -> int:
127 """Count the number of pixels in an image with a given mask bit set.
129 Parameters
130 ----------
131 maskedImage : `lsst.afw.image.MaskedImage`
132 The masked image,
133 maskPlane : `str`
134 The name of the bitmask.
136 Returns
137 -------
138 count : `int``
139 The number of pixels in with the selected mask bit
140 """
141 bit = maskedImage.mask.getPlaneBitMask(maskPlane)
142 return len(np.where(np.bitwise_and(maskedImage.mask.array, bit))[0])
145def quickSmooth(data: npt.NDArray[np.float64], sigma: float = 2) -> npt.NDArray[np.float64]:
146 """Perform a quick smoothing of the image.
148 Not to be used for scientific purposes, but improves the stretch and
149 visual rendering of low SNR against the sky background in cutouts.
151 Parameters
152 ----------
153 data : `np.array`
154 The image data to smooth
155 sigma : `float`, optional
156 The size of the smoothing kernel.
158 Returns
159 -------
160 smoothData : `np.array`
161 The smoothed data
162 """
163 kernel = [sigma, sigma]
164 smoothData = gaussian_filter(data, kernel, mode="constant")
165 return smoothData
168def argMax2d(array: npt.NDArray[np.float64]) -> tuple[tuple[float, float], bool, list[tuple[float, float]]]:
169 """Get the index of the max value of an array and whether it's unique.
171 If its not unique, returns a list of the other locations containing the
172 maximum value, e.g. returns
174 (12, 34), False, [(56,78), (910, 1112)]
176 Parameters
177 ----------
178 array : `np.array`
179 The data
181 Returns
182 -------
183 maxLocation : `tuple`
184 The coords of the first instance of the max value
185 unique : `bool`
186 Whether it's the only location
187 otherLocations : `list` of `tuple`
188 List of the other max values' locations, empty if False
189 """
190 uniqueMaximum = False
191 maxCoords = np.where(array == np.max(array))
192 # list of coords as tuples
193 listMaxCoords = [(float(coord[0]), float(coord[1])) for coord in zip(*maxCoords)]
194 if len(listMaxCoords) == 1: # single unambiguous value
195 uniqueMaximum = True
197 return listMaxCoords[0], uniqueMaximum, listMaxCoords[1:]
200def dayObsIntToString(dayObs: int) -> str:
201 """Convert an integer dayObs to a dash-delimited string.
203 e.g. convert the hard to read 20210101 to 2021-01-01
205 Parameters
206 ----------
207 dayObs : `int`
208 The dayObs.
210 Returns
211 -------
212 dayObs : `str`
213 The dayObs as a string.
214 """
215 assert isinstance(dayObs, int)
216 dStr = str(dayObs)
217 assert len(dStr) == 8
218 return "-".join([dStr[0:4], dStr[4:6], dStr[6:8]])
221def dayObsSeqNumToVisitId(dayObs: int, seqNum: int) -> int:
222 """Get the visit id for a given dayObs/seqNum.
224 Parameters
225 ----------
226 dayObs : `int`
227 The dayObs.
228 seqNum : `int`
229 The seqNum.
231 Returns
232 -------
233 visitId : `int`
234 The visitId.
236 Notes
237 -----
238 TODO: Remove this horrible hack once DM-30948 makes this possible
239 programatically/via the butler.
240 """
241 if dayObs < 19700101 or dayObs > 35000101:
242 raise ValueError(f"dayObs value {dayObs} outside plausible range")
243 return int(f"{dayObs}{seqNum:05}")
246def getImageStats(exp: afwImage.Exposure) -> pipeBase.Struct:
247 """Calculate a grab-bag of stats for an image. Must remain fast.
249 Parameters
250 ----------
251 exp : `lsst.afw.image.Exposure`
252 The input exposure.
254 Returns
255 -------
256 stats : `lsst.pipe.base.Struct`
257 A container with attributes containing measurements and statistics
258 for the image.
259 """
260 result = pipeBase.Struct()
262 vi = exp.visitInfo
263 expTime = vi.exposureTime
264 md = exp.getMetadata()
266 obj = vi.object
267 mjd = vi.getDate().get()
268 result.object = obj
269 result.mjd = mjd
271 fullFilterString = exp.filter.physicalLabel
272 filt = fullFilterString.split(FILTER_DELIMITER)[0]
273 grating = fullFilterString.split(FILTER_DELIMITER)[1]
275 airmass = vi.getBoresightAirmass()
276 rotangle = vi.getBoresightRotAngle().asDegrees()
278 azAlt = vi.getBoresightAzAlt()
279 az = azAlt[0].asDegrees()
280 el = azAlt[1].asDegrees()
282 result.expTime = expTime
283 result.filter = filt
284 result.grating = grating
285 result.airmass = airmass
286 result.rotangle = rotangle
287 result.az = az
288 result.el = el
289 result.focus = md.get("FOCUSZ")
291 data = exp.image.array
292 result.maxValue = np.max(data)
294 peak, uniquePeak, otherPeaks = argMax2d(data)
295 result.maxPixelLocation = peak
296 result.multipleMaxPixels = uniquePeak
298 result.nBadPixels = countPixels(exp.maskedImage, "BAD")
299 result.nSatPixels = countPixels(exp.maskedImage, "SAT")
300 result.percentile99 = np.percentile(data, 99)
301 result.percentile9999 = np.percentile(data, 99.99)
303 sctrl = afwMath.StatisticsControl()
304 sctrl.setNumSigmaClip(5)
305 sctrl.setNumIter(2)
306 statTypes = afwMath.MEANCLIP | afwMath.STDEVCLIP
307 stats = afwMath.makeStatistics(exp.maskedImage, statTypes, sctrl)
308 std, stderr = stats.getResult(afwMath.STDEVCLIP)
309 mean, meanerr = stats.getResult(afwMath.MEANCLIP)
311 result.clippedMean = mean
312 result.clippedStddev = std
314 return result
317def detectObjectsInExp(
318 exp: afwImage.Exposure, nSigma: float = 10, nPixMin: int = 10, grow: int = 0
319) -> afwDetect.FootprintSet:
320 """Quick and dirty object detection for an exposure.
322 Return the footPrintSet for the objects in a preferably-postISR exposure.
324 Parameters
325 ----------
326 exp : `lsst.afw.image.Exposure`
327 The exposure to detect objects in.
328 nSigma : `float`
329 The number of sigma for detection.
330 nPixMin : `int`
331 The minimum number of pixels in an object for detection.
332 grow : `int`
333 The number of pixels to grow the footprint by after detection.
335 Returns
336 -------
337 footPrintSet : `lsst.afw.detection.FootprintSet`
338 The set of footprints in the image.
339 """
340 median = np.nanmedian(exp.image.array)
341 exp.image -= median
343 threshold = afwDetect.Threshold(nSigma, afwDetect.Threshold.STDEV)
344 footPrintSet = afwDetect.FootprintSet(exp.getMaskedImage(), threshold, "DETECTED", nPixMin)
345 if grow > 0:
346 isotropic = True
347 footPrintSet = afwDetect.FootprintSet(footPrintSet, grow, isotropic)
349 exp.image += median # add back in to leave background unchanged
350 return footPrintSet
353def fluxesFromFootprints(
354 footprints: afwDetect.FootprintSet | afwDetect.Footprint | Iterable[afwDetect.Footprint],
355 parentImage: afwImage.Image,
356 subtractImageMedian: bool = False,
357) -> npt.NDArray[np.float64]:
358 """Calculate the flux from a set of footprints, given the parent image,
359 optionally subtracting the whole-image median from each pixel as a very
360 rough background subtraction.
362 Parameters
363 ----------
364 footprints : `lsst.afw.detection.FootprintSet` or
365 `lsst.afw.detection.Footprint` or
366 `iterable` of `lsst.afw.detection.Footprint`
367 The footprints to measure.
368 parentImage : `lsst.afw.image.Image`
369 The parent image.
370 subtractImageMedian : `bool`, optional
371 Subtract a whole-image median from each pixel in the footprint when
372 summing as a very crude background subtraction. Does not change the
373 original image.
375 Returns
376 -------
377 fluxes : `list` of `float`
378 The fluxes for each footprint.
380 Raises
381 ------
382 TypeError : raise for unsupported types.
383 """
384 median = 0
385 if subtractImageMedian:
386 median = np.nanmedian(parentImage.array)
388 # poor person's single dispatch
389 badTypeMsg = (
390 "This function works with FootprintSets, single Footprints, and iterables of Footprints. "
391 f"Got {type(footprints)}: {footprints}"
392 )
393 if isinstance(footprints, FootprintSet):
394 footprints = footprints.getFootprints()
395 elif isinstance(footprints, Sequence):
396 if not isinstance(footprints[0], Footprint):
397 raise TypeError(badTypeMsg)
398 elif isinstance(footprints, Footprint):
399 footprints = [footprints]
400 else:
401 raise TypeError(badTypeMsg)
403 return np.array([fluxFromFootprint(fp, parentImage, backgroundValue=median) for fp in footprints])
406def fluxFromFootprint(
407 footprint: afwDetection.Footprint, parentImage: afwImage.Image, backgroundValue: float = 0
408) -> float:
409 """Calculate the flux from a footprint, given the parent image, optionally
410 subtracting a single value from each pixel as a very rough background
411 subtraction, e.g. the image median.
413 Parameters
414 ----------
415 footprint : `lsst.afw.detection.Footprint`
416 The footprint to measure.
417 parentImage : `lsst.afw.image.Image`
418 Image containing the footprint.
419 backgroundValue : `bool`, optional
420 The value to subtract from each pixel in the footprint when summing
421 as a very crude background subtraction. Does not change the original
422 image.
424 Returns
425 -------
426 flux : `float`
427 The flux in the footprint
428 """
429 if backgroundValue: # only do the subtraction if non-zero for speed
430 xy0 = parentImage.getBBox().getMin()
431 return footprint.computeFluxFromArray(parentImage.array - backgroundValue, xy0)
432 return footprint.computeFluxFromImage(parentImage)
435def humanNameForCelestialObject(objName: str) -> list[str]:
436 """Returns a list of all human names for obj, or [] if none are found.
438 Parameters
439 ----------
440 objName : `str`
441 The/a name of the object.
443 Returns
444 -------
445 names : `list` of `str`
446 The names found for the object
447 """
448 from astroquery.simbad import Simbad
450 results = []
451 try:
452 simbadResult = Simbad.query_objectids(objName)
453 for row in simbadResult:
454 if row["ID"].startswith("NAME"):
455 results.append(row["ID"].replace("NAME ", ""))
456 return results
457 except Exception:
458 return [] # same behavior as for found but un-named objects
461def _getAltAzZenithsFromSeqNum(
462 butler: dafButler.Butler, dayObs: int, seqNumList: list[int]
463) -> tuple[list[float], list[float], list[float]]:
464 """Get the alt, az and zenith angle for the seqNums of a given dayObs.
466 Parameters
467 ----------
468 butler : `lsst.daf.butler.Butler`
469 The butler to query.
470 dayObs : `int`
471 The dayObs.
472 seqNumList : `list` of `int`
473 The seqNums for which to return the alt, az and zenith
475 Returns
476 -------
477 azimuths : `list` of `float`
478 List of the azimuths for each seqNum
479 elevations : `list` of `float`
480 List of the elevations for each seqNum
481 zeniths : `list` of `float`
482 List of the zenith angles for each seqNum
483 """
484 azimuths, elevations, zeniths = [], [], []
485 for seqNum in seqNumList:
486 md = butler.get("raw.metadata", day_obs=dayObs, seq_num=seqNum, detector=0)
487 obsInfo = ObservationInfo(md)
488 alt = obsInfo.altaz_begin.alt.value
489 az = obsInfo.altaz_begin.az.value
490 elevations.append(alt)
491 zeniths.append(90 - alt)
492 azimuths.append(az)
493 return azimuths, elevations, zeniths
496def getFocusFromHeader(exp: afwImage.Exposure) -> float | None:
497 """Get the raw focus value from the header.
499 Parameters
500 ----------
501 exp : `lsst.afw.image.exposure`
502 The exposure.
504 Returns
505 -------
506 focus : `float` or `None`
507 The focus value if found, else ``None``.
508 """
509 md = exp.getMetadata()
510 if "FOCUSZ" in md:
511 return md["FOCUSZ"]
512 return None
515def checkStackSetup() -> None:
516 """Check which weekly tag is being used and which local packages are setup.
518 Designed primarily for use in notbooks/observing, this prints the weekly
519 tag(s) are setup for lsst_distrib, and lists any locally setup packages and
520 the path to each.
522 Notes
523 -----
524 Uses print() instead of logger messages as this should simply print them
525 without being vulnerable to any log messages potentially being diverted.
526 """
527 packages = packageUtils.getEnvironmentPackages(include_all=True)
529 lsstDistribHashAndTags = packages["lsst_distrib"] # looks something like 'g4eae7cb9+1418867f (w_2022_13)'
530 lsstDistribTags = lsstDistribHashAndTags.split()[1]
531 if len(lsstDistribTags.split()) == 1:
532 tag = lsstDistribTags.replace("(", "")
533 tag = tag.replace(")", "")
534 print(f"You are running {tag} of lsst_distrib")
535 else: # multiple weekly tags found for lsst_distrib!
536 print(f"The version of lsst_distrib you have is compatible with: {lsstDistribTags}")
538 localPackages = []
539 localPaths = []
540 for package, tags in packages.items():
541 if tags.startswith("LOCAL:"):
542 path = tags.split("LOCAL:")[1]
543 path = path.split("@")[0] # don't need the git SHA etc
544 localPaths.append(path)
545 localPackages.append(package)
547 if localPackages:
548 print("\nLocally setup packages:")
549 print("-----------------------")
550 maxLen = max(len(package) for package in localPackages)
551 for package, path in zip(localPackages, localPaths):
552 print(f"{package:<{maxLen}s} at {path}")
553 else:
554 print("\nNo locally setup packages (using a vanilla stack)")
557def setupLogging(longlog: bool = False) -> None:
558 """Setup logging in the same way as one would get from pipetask run.
560 Code that isn't run through the butler CLI defaults to WARNING level
561 messages and no logger names. This sets the behaviour to follow whatever
562 the pipeline default is, currently
563 <logger_name> <level>: <message> e.g.
564 lsst.isr INFO: Masking defects.
565 """
566 CliLog.initLog(longlog=longlog)
569def getCurrentDayObs_datetime() -> datetime.date:
570 """Get the current day_obs - the observatory rolls the date over at UTC-12
572 Returned as datetime.date(2022, 4, 28)
573 """
574 utc = gettz("UTC")
575 nowUtc = datetime.datetime.now().astimezone(utc)
576 offset = datetime.timedelta(hours=-12)
577 dayObs = (nowUtc + offset).date()
578 return dayObs
581def getCurrentDayObs_int() -> int:
582 """Return the current dayObs as an int in the form 20220428"""
583 return int(getCurrentDayObs_datetime().strftime("%Y%m%d"))
586def getCurrentDayObs_humanStr() -> str:
587 """Return the current dayObs as a string in the form '2022-04-28'"""
588 return dayObsIntToString(getCurrentDayObs_int())
591def getExpRecordAge(expRecord: dafButler.DimensionRecord) -> float:
592 """Get the time, in seconds, since the end of exposure.
594 Parameters
595 ----------
596 expRecord : `lsst.daf.butler.DimensionRecord`
597 The exposure record.
599 Returns
600 -------
601 age : `float`
602 The age of the exposure, in seconds.
603 """
604 return -1 * (expRecord.timespan.end - Time.now()).sec
607def getSite() -> str:
608 """Returns where the code is running.
610 Returns
611 -------
612 location : `str`
613 One of:
614 ['tucson', 'summit', 'base', 'staff-rsp', 'rubin-devl', 'jenkins',
615 'usdf-k8s']
616 If the location cannot be determined "UNKOWN" is returned.
618 """
619 # All nublado instances guarantee that EXTERNAL_URL is set and uniquely
620 # identifies it.
621 location = os.getenv("EXTERNAL_INSTANCE_URL", "")
622 if location == "https://tucson-teststand.lsst.codes": 622 ↛ 623line 622 didn't jump to line 623 because the condition on line 622 was never true
623 return "tucson"
624 elif location == "https://summit-lsp.lsst.codes": 624 ↛ 625line 624 didn't jump to line 625 because the condition on line 624 was never true
625 return "summit"
626 elif location == "https://base-lsp.lsst.codes": 626 ↛ 627line 626 didn't jump to line 627 because the condition on line 626 was never true
627 return "base"
628 elif location == "https://usdf-rsp.slac.stanford.edu": 628 ↛ 629line 628 didn't jump to line 629 because the condition on line 628 was never true
629 return "staff-rsp"
630 elif location == "https://usdf-rsp-dev.slac.stanford.edu": 630 ↛ 631line 630 didn't jump to line 631 because the condition on line 630 was never true
631 return "staff-rsp" # we don't care this is the dev RSP, it's basically the same env wrt paths etc
633 # if no EXTERNAL_URL, try HOSTNAME to see if we're on the dev nodes
634 # it is expected that this will be extensible to SLAC
635 hostname = os.getenv("HOSTNAME", "")
636 interactiveNodes = ("sdfrome", "sdfiana")
637 if hostname.startswith(interactiveNodes): 637 ↛ 638line 637 didn't jump to line 638 because the condition on line 637 was never true
638 return "rubin-devl"
639 elif hostname == "htcondor.ls.lsst.org": 639 ↛ 640line 639 didn't jump to line 640 because the condition on line 639 was never true
640 return "base"
641 elif hostname == "htcondor.cp.lsst.org": 641 ↛ 642line 641 didn't jump to line 642 because the condition on line 641 was never true
642 return "summit"
644 jenkinsHome = os.getenv("JENKINS_HOME", "")
645 if jenkinsHome != "": 645 ↛ 649line 645 didn't jump to line 649 because the condition on line 645 was always true
646 return "jenkins"
648 # we're probably inside a k8s pod doing rapid analysis work at this point
649 location = os.getenv("RAPID_ANALYSIS_LOCATION", "")
650 if location == "TTS":
651 return "tucson"
652 if location == "BTS":
653 return "base"
654 if location == "SUMMIT":
655 return "summit"
656 if location == "USDF":
657 return "usdf-k8s"
659 # we have failed
660 return "UNKOWN"
663def getAltAzFromSkyPosition(
664 skyPos: geom.SpherePoint,
665 visitInfo: afwImage.VisitInfo,
666 doCorrectRefraction: bool = False,
667 wavelength: float = 500.0,
668 pressureOverride: float | None = None,
669 temperatureOverride: float | None = None,
670 relativeHumidityOverride: float | None = None,
671) -> tuple[geom.Angle, geom.Angle]:
672 """Get the alt/az from the position on the sky and the time and location
673 of the observation.
675 The temperature, pressure and relative humidity are taken from the
676 visitInfo by default, but can be individually overridden as needed. It
677 should be noted that the visitInfo never contains a nominal wavelength, and
678 so this takes a default value of 500nm.
680 Parameters
681 ----------
682 skyPos : `lsst.geom.SpherePoint`
683 The position on the sky.
684 visitInfo : `lsst.afw.image.VisitInfo`
685 The visit info containing the time of the observation.
686 doCorrectRefraction : `bool`, optional
687 Correct for the atmospheric refraction?
688 wavelength : `float`, optional
689 The nominal wavelength in nanometers (e.g. 500.0), as a float.
690 pressureOverride : `float`, optional
691 The pressure, in bars (e.g. 0.770), to override the value supplied in
692 the visitInfo, as a float.
693 temperatureOverride : `float`, optional
694 The temperature, in Celsius (e.g. 10.0), to override the value supplied
695 in the visitInfo, as a float.
696 relativeHumidityOverride : `float`, optional
697 The relativeHumidity in the range 0..1 (i.e. not as a percentage), to
698 override the value supplied in the visitInfo, as a float.
700 Returns
701 -------
702 alt : `lsst.geom.Angle`
703 The altitude.
704 az : `lsst.geom.Angle`
705 The azimuth.
706 """
707 skyLocation = SkyCoord(skyPos.getRa().asRadians(), skyPos.getDec().asRadians(), unit=u.rad)
708 long = visitInfo.observatory.getLongitude()
709 lat = visitInfo.observatory.getLatitude()
710 ele = visitInfo.observatory.getElevation()
711 earthLocation = EarthLocation.from_geodetic(long.asDegrees(), lat.asDegrees(), ele)
713 refractionKwargs = {}
714 if doCorrectRefraction:
715 # wavelength is never supplied in the visitInfo so always take this
716 wavelength = wavelength * u.nm
718 if pressureOverride:
719 pressure = pressureOverride
720 else:
721 pressure = visitInfo.weather.getAirPressure()
722 # ObservationInfos (which are the "source of truth" use pascals) so
723 # convert from pascals to bars
724 pressure /= 100000.0
725 pressure = pressure * u.bar
727 if temperatureOverride:
728 temperature = temperatureOverride
729 else:
730 temperature = visitInfo.weather.getAirTemperature()
731 temperature = temperature * u.deg_C
733 if relativeHumidityOverride:
734 relativeHumidity = relativeHumidityOverride
735 else:
736 relativeHumidity = visitInfo.weather.getHumidity() / 100.0 # this is in percent
737 relativeHumidity = relativeHumidity
739 refractionKwargs = dict(
740 pressure=pressure, temperature=temperature, relative_humidity=relativeHumidity, obswl=wavelength
741 )
743 obsTime = visitInfo.date.toAstropy()
744 altAz = AltAz(obstime=obsTime, location=earthLocation, **refractionKwargs)
746 obsAltAz = skyLocation.transform_to(altAz)
747 alt = geom.Angle(obsAltAz.alt.degree, geom.degrees)
748 az = geom.Angle(obsAltAz.az.degree, geom.degrees)
750 return alt, az
753def getExpPositionOffset(
754 exp1: afwImage.Exposure,
755 exp2: afwImage.Exposure,
756 useWcs: bool = True,
757 allowDifferentPlateScales: bool = False,
758) -> pipeBase.Struct:
759 """Get the change in sky position between two exposures.
761 Given two exposures, calculate the offset on the sky between the images.
762 If useWcs then use the (fitted or unfitted) skyOrigin from their WCSs, and
763 calculate the alt/az from the observation times, otherwise use the nominal
764 values in the exposures' visitInfos. Note that if using the visitInfo
765 values that for a given pointing the ra/dec will be ~identical, regardless
766 of whether astrometric fitting has been performed.
768 Values are given as exp1-exp2.
770 Parameters
771 ----------
772 exp1 : `lsst.afw.image.Exposure`
773 The first exposure.
774 exp2 : `lsst.afw.image.Exposure`
775 The second exposure.
776 useWcs : `bool`
777 Use the WCS for the ra/dec and alt/az if True, else use the nominal/
778 boresight values from the exposures' visitInfos.
779 allowDifferentPlateScales : `bool`, optional
780 Use to disable checking that plate scales are the same. Generally,
781 differing plate scales would indicate an error, but where blind-solving
782 has been undertaken during commissioning plate scales can be different
783 enough to warrant setting this to ``True``.
785 Returns
786 -------
787 offsets : `lsst.pipe.base.Struct`
788 A struct containing the offsets:
789 ``deltaRa``
790 The diference in ra (`lsst.geom.Angle`)
791 ``deltaDec``
792 The diference in dec (`lsst.geom.Angle`)
793 ``deltaAlt``
794 The diference in alt (`lsst.geom.Angle`)
795 ``deltaAz``
796 The diference in az (`lsst.geom.Angle`)
797 ``deltaPixels``
798 The diference in pixels (`float`)
799 """
801 wcs1 = exp1.getWcs()
802 wcs2 = exp2.getWcs()
803 pixScaleArcSec = wcs1.getPixelScale().asArcseconds()
804 if not allowDifferentPlateScales:
805 assert np.isclose(
806 pixScaleArcSec, wcs2.getPixelScale().asArcseconds()
807 ), "Pixel scales in the exposures differ."
809 if useWcs:
810 p1 = wcs1.getSkyOrigin()
811 p2 = wcs2.getSkyOrigin()
812 alt1, az1 = getAltAzFromSkyPosition(p1, exp1.getInfo().getVisitInfo())
813 alt2, az2 = getAltAzFromSkyPosition(p2, exp2.getInfo().getVisitInfo())
814 ra1 = p1[0]
815 ra2 = p2[0]
816 dec1 = p1[1]
817 dec2 = p2[1]
818 else:
819 az1 = exp1.visitInfo.boresightAzAlt[0]
820 az2 = exp2.visitInfo.boresightAzAlt[0]
821 alt1 = exp1.visitInfo.boresightAzAlt[1]
822 alt2 = exp2.visitInfo.boresightAzAlt[1]
824 ra1 = exp1.visitInfo.boresightRaDec[0]
825 ra2 = exp2.visitInfo.boresightRaDec[0]
826 dec1 = exp1.visitInfo.boresightRaDec[1]
827 dec2 = exp2.visitInfo.boresightRaDec[1]
829 p1 = exp1.visitInfo.boresightRaDec
830 p2 = exp2.visitInfo.boresightRaDec
832 angular_offset = p1.separation(p2).asArcseconds()
833 deltaPixels = angular_offset / pixScaleArcSec
835 ret = pipeBase.Struct(
836 deltaRa=(ra1 - ra2).wrapNear(geom.Angle(0.0)),
837 deltaDec=dec1 - dec2,
838 deltaAlt=alt1 - alt2,
839 deltaAz=(az1 - az2).wrapNear(geom.Angle(0.0)),
840 deltaPixels=deltaPixels,
841 )
843 return ret
846def starTrackerFileToExposure(filename: str, logger: logging.Logger | None = None) -> afwImage.Exposure:
847 """Read the exposure from the file and set the wcs from the header.
849 Parameters
850 ----------
851 filename : `str`
852 The full path to the file.
853 logger : `logging.Logger`, optional
854 The logger to use for errors, created if not supplied.
856 Returns
857 -------
858 exp : `lsst.afw.image.Exposure`
859 The exposure.
860 """
861 if not logger:
862 logger = logging.getLogger(__name__)
863 exp = afwImage.ExposureF(filename)
864 try:
865 wcs = genericCameraHeaderToWcs(exp)
866 exp.setWcs(wcs)
867 except Exception as e:
868 logger.warning(f"Failed to set wcs from header: {e}")
870 # for some reason the date isn't being set correctly
871 # DATE-OBS is present in the original header, but it's being
872 # stripped out and somehow not set (plus it doesn't give the midpoint
873 # of the exposure), so set it manually from the midpoint here
874 try:
875 newArgs = {} # dict to unpack into visitInfo.copyWith - fill it with whatever needs to be replaced
876 md = exp.getMetadata()
878 begin = datetime.datetime.fromisoformat(md["DATE-BEG"])
879 end = datetime.datetime.fromisoformat(md["DATE-END"])
880 duration = end - begin
881 mid = begin + duration / 2
882 newTime = dafBase.DateTime(mid.isoformat(), dafBase.DateTime.Timescale.TAI)
883 newArgs["date"] = newTime
885 # AIRPRESS is being set as PRESSURE so afw doesn't pick it up
886 # once we're using the butler for data we will just set it to take
887 # PRESSURE in the translator instead of this
888 weather = exp.visitInfo.getWeather()
889 oldPressure = weather.getAirPressure()
890 if not np.isfinite(oldPressure):
891 pressure = md.get("PRESSURE")
892 if pressure is not None:
893 logger.info("Patching the weather info using the PRESSURE header keyword")
894 newWeather = Weather(weather.getAirTemperature(), pressure, weather.getHumidity())
895 newArgs["weather"] = newWeather
897 if newArgs:
898 newVi = exp.visitInfo.copyWith(**newArgs)
899 exp.info.setVisitInfo(newVi)
900 except Exception as e:
901 logger.warning(f"Failed to set date from header: {e}")
903 return exp
906def obsInfoToDict(obsInfo: ObservationInfo) -> dict:
907 """Convert an ObservationInfo to a dict.
909 Parameters
910 ----------
911 obsInfo : `astro_metadata_translator.ObservationInfo`
912 The ObservationInfo to convert.
914 Returns
915 -------
916 obsInfoDict : `dict`
917 The ObservationInfo as a dict.
918 """
919 return {prop: getattr(obsInfo, prop) for prop in obsInfo.all_properties.keys()}
922def getFieldNameAndTileNumber(
923 field: str, warn: bool = True, logger: logging.Logger | None = None
924) -> tuple[str, int | None]:
925 """Get the tile name and number of an observed field.
927 It is assumed to always be appended, with an underscore, to the rest of the
928 field name. Returns the name and number as a tuple, or the name unchanged
929 if no tile number is found.
931 Parameters
932 ----------
933 field : `str`
934 The name of the field
936 Returns
937 -------
938 fieldName : `str`
939 The name of the field without the trailing tile number, if present.
940 tileNum : `int`
941 The number of the tile, as an integer, or ``None`` if not found.
942 """
943 if warn and not logger:
944 logger = logging.getLogger("lsst.summit.utils.utils.getFieldNameAndTileNumber")
946 if "_" not in field:
947 if warn and logger is not None:
948 logger.warning(
949 f"Field {field} does not contain an underscore," " so cannot determine the tile number."
950 )
951 return field, None
953 try:
954 fieldParts = field.split("_")
955 fieldNum = int(fieldParts[-1])
956 except ValueError:
957 if warn and logger is not None:
958 logger.warning(
959 f"Field {field} does not contain only an integer after the final underscore"
960 " so cannot determine the tile number."
961 )
962 return field, None
964 return "_".join(fieldParts[:-1]), fieldNum
967def getAirmassSeeingCorrection(airmass: float) -> float:
968 """Get the correction factor for seeing due to airmass.
970 Parameters
971 ----------
972 airmass : `float`
973 The airmass, greater than or equal to 1.
975 Returns
976 -------
977 correctionFactor : `float`
978 The correction factor to apply to the seeing.
980 Raises
981 ------
982 ValueError raised for unphysical airmasses.
983 """
984 if airmass < 1:
985 raise ValueError(f"Invalid airmass: {airmass}")
986 return airmass ** (-0.6)
989@deprecated(
990 reason=". Will be removed after v28.0.",
991 version="v27.0",
992 category=FutureWarning,
993)
994def getFilterSeeingCorrection(filterName: str) -> float:
995 """Get the correction factor for seeing due to a filter.
997 Parameters
998 ----------
999 filterName : `str`
1000 The name of the filter, e.g. 'SDSSg_65mm'.
1002 Returns
1003 -------
1004 correctionFactor : `float`
1005 The correction factor to apply to the seeing.
1007 Raises
1008 ------
1009 ValueError raised for unknown filters.
1010 """
1011 return getBandpassSeeingCorrection(filterName)
1014def getBandpassSeeingCorrection(filterName: str) -> float:
1015 """Get the correction factor for seeing due to a filter.
1017 Parameters
1018 ----------
1019 filterName : `str`
1020 The name of the filter, e.g. 'SDSSg_65mm'.
1022 Returns
1023 -------
1024 correctionFactor : `float`
1025 The correction factor to apply to the seeing.
1027 Raises
1028 ------
1029 ValueError raised for unknown filters.
1030 """
1031 match filterName:
1032 case "SDSSg_65mm": # LATISS
1033 return (474.41 / 500.0) ** 0.2
1034 case "SDSSr_65mm": # LATISS
1035 return (628.47 / 500.0) ** 0.2
1036 case "SDSSi_65mm": # LATISS
1037 return (769.51 / 500.0) ** 0.2
1038 case "SDSSz_65mm": # LATISS
1039 return (871.45 / 500.0) ** 0.2
1040 case "SDSSy_65mm": # LATISS
1041 return (986.8 / 500.0) ** 0.2
1042 case "u_02": # ComCam
1043 return (370.697 / 500.0) ** 0.2
1044 case "g_01": # ComCam
1045 return (476.359 / 500.0) ** 0.2
1046 case "r_03": # ComCam
1047 return (619.383 / 500.0) ** 0.2
1048 case "i_06": # ComCam
1049 return (754.502 / 500.0) ** 0.2
1050 case "z_03": # ComCam
1051 return (866.976 / 500.0) ** 0.2
1052 case "y_04": # ComCam
1053 return (972.713 / 500.0) ** 0.2
1054 case "u_24": # LSSTCam
1055 return (365.1397 / 500.0) ** 0.2
1056 case "g_6": # LSSTCam
1057 return (475.6104 / 500.0) ** 0.2
1058 case "r_57": # LSSTCam
1059 return (618.6055 / 500.0) ** 0.2
1060 case "i_39": # LSSTCam
1061 return (753.5490 / 500.0) ** 0.2
1062 case "z_20": # LSSTCam
1063 return (868.0935 / 500.0) ** 0.2
1064 case "y_10": # LSSTCam
1065 return (972.5979 / 500.0) ** 0.2
1066 case _:
1067 raise ValueError(f"Unknown filter name: {filterName}")
1070def getCdf(data: np.ndarray, scale: int, nBinsMax: int = 300_000) -> tuple[np.ndarray | float, float, float]:
1071 """Return an approximate cumulative distribution function scaled to
1072 the [0, scale] range.
1074 If the input data is all nan, then the output cdf will be nan as well as
1075 the min and max values.
1077 Parameters
1078 ----------
1079 data : `np.array`
1080 The input data.
1081 scale : `int`
1082 The scaling range of the output.
1083 nBinsMax : `int`, optional
1084 Maximum number of bins to use.
1086 Returns
1087 -------
1088 cdf : `np.array` of `int`
1089 A monotonically increasing sequence that represents a scaled
1090 cumulative distribution function, starting with the value at
1091 minVal, then at (minVal + 1), and so on.
1092 minVal : `float`
1093 An integer smaller than the minimum value in the input data.
1094 maxVal : `float`
1095 An integer larger than the maximum value in the input data.
1096 """
1097 flatData = data.ravel()
1098 size = flatData.size - np.count_nonzero(np.isnan(flatData))
1100 minVal = np.floor(np.nanmin(flatData))
1101 maxVal = np.ceil(np.nanmax(flatData)) + 1.0
1103 if np.isnan(minVal) or np.isnan(maxVal):
1104 # if either the min or max are nan, then the data is all nan as we're
1105 # using nanmin and nanmax. Given this, we can't calculate a cdf, so
1106 # return nans for all values
1107 return np.nan, np.nan, np.nan
1109 nBins = np.clip(int(maxVal) - int(minVal), 1, nBinsMax)
1111 hist, binEdges = np.histogram(flatData, bins=nBins, range=(int(minVal), int(maxVal)))
1113 cdf = (scale * np.cumsum(hist) / size).astype(np.int64)
1114 return cdf, minVal, maxVal
1117def getQuantiles(data: npt.NDArray[np.float64], nColors: int) -> npt.NDArray[np.float64]:
1118 """Get a set of boundaries that equally distribute data into
1119 nColors intervals. The output can be used to make a colormap of nColors
1120 colors.
1122 This is equivalent to using the numpy function:
1123 np.nanquantile(data, np.linspace(0, 1, nColors + 1))
1124 but with a coarser precision, yet sufficient for our use case. This
1125 implementation gives a significant speed-up. In the case of large
1126 ranges, np.nanquantile is used because it is more memory efficient.
1128 If all elements of ``data`` are nan then the output ``boundaries`` will
1129 also all be ``nan`` to keep the interface consistent.
1131 Parameters
1132 ----------
1133 data : `np.array`
1134 The input image data.
1135 nColors : `int`
1136 The number of intervals to distribute data into.
1138 Returns
1139 -------
1140 boundaries: `list` of `float`
1141 A monotonically increasing sequence of size (nColors + 1). These are
1142 the edges of nColors intervals.
1143 """
1144 dataRange = np.nanmax(data) - np.nanmin(data)
1145 if dataRange > 300_000:
1146 # Use slower but memory efficient nanquantile
1147 logger = logging.getLogger(__name__)
1148 logger.warning(f"Data range is very large ({dataRange}); using slower quantile code.")
1149 boundaries = np.nanquantile(data, np.linspace(0, 1, nColors + 1))
1150 else:
1151 cdf, minVal, maxVal = getCdf(data, nColors)
1152 if np.isnan(minVal): # cdf calculation has failed because all data is nan
1153 return np.asarray([np.nan for _ in range(nColors)])
1154 assert isinstance(cdf, np.ndarray), "cdf is not an np.ndarray"
1155 scale = (maxVal - minVal) / len(cdf)
1157 boundaries = np.asarray([np.argmax(cdf >= i) * scale + minVal for i in range(nColors)] + [maxVal])
1159 return boundaries
1162def digitizeData(data: npt.NDArray[np.float64], nColors: int = 256) -> npt.NDArray[np.integer]:
1163 """
1164 Scale data into nColors using its cumulative distribution function.
1166 Parameters
1167 ----------
1168 data : `np.array`
1169 The input image data.
1170 nColors : `int`
1171 The number of intervals to distribute data into.
1173 Returns
1174 -------
1175 data: `np.array` of `int`
1176 Scaled data in the [0, nColors - 1] range.
1177 """
1178 cdf, minVal, maxVal = getCdf(data, nColors - 1)
1179 assert isinstance(cdf, np.ndarray), "cdf is not an np.ndarray"
1180 scale = (maxVal - minVal) / len(cdf)
1181 bins = np.floor((data * scale - minVal)).astype(np.int64)
1182 return cdf[bins]
1185def getBboxAround(centroid: geom.Point, boxSize: int, exp: afwImage.Exposure) -> geom.Box2I:
1186 """Get a bbox centered on a point, clipped to the exposure.
1188 If the bbox would extend beyond the bounds of the exposure it is clipped to
1189 the exposure, resulting in a non-square bbox.
1191 Parameters
1192 ----------
1193 centroid : `lsst.geom.Point`
1194 The source centroid.
1195 boxsize : `int`
1196 The size of the box to centre around the centroid.
1197 exp : `lsst.afw.image.Exposure`
1198 The exposure, so the bbox can be clipped to not overrun the bounds.
1200 Returns
1201 -------
1202 bbox : `lsst.geom.Box2I`
1203 The bounding box, centered on the centroid unless clipping to the
1204 exposure causes it to be non-square.
1205 """
1206 bbox = geom.BoxI().makeCenteredBox(centroid, geom.Extent2I(boxSize, boxSize))
1207 bbox = bbox.clippedTo(exp.getBBox())
1208 return bbox
1211def bboxToMatplotlibRectanle(bbox: geom.Box2I | geom.Box2D) -> matplotlib.patches.Rectangle:
1212 """Convert a bbox to a matplotlib Rectangle for plotting.
1214 Parameters
1215 ----------
1216 results : `lsst.geom.Box2I` or `lsst.geom.Box2D`
1217 The bbox to convert.
1219 Returns
1220 -------
1221 rectangle : `matplotlib.patches.Rectangle`
1222 The rectangle.
1223 """
1224 ll = bbox.minX, bbox.minY
1225 width, height = bbox.getDimensions()
1226 return Rectangle(ll, width, height)
1229def computeExposureId(instrument: str, controller: str, dayObs: int, seqNum: int) -> int:
1230 instrument = instrument.lower()
1231 if instrument == "latiss":
1232 return lsst.obs.lsst.translators.LatissTranslator.compute_exposure_id(dayObs, seqNum, controller)
1233 elif instrument == "lsstcomcam":
1234 return lsst.obs.lsst.translators.LsstComCamTranslator.compute_exposure_id(dayObs, seqNum, controller)
1235 elif instrument == "lsstcomcamsim":
1236 return lsst.obs.lsst.translators.LsstComCamSimTranslator.compute_exposure_id(
1237 dayObs, seqNum, controller
1238 )
1239 elif instrument == "lsstcam":
1240 return lsst.obs.lsst.translators.LsstCamTranslator.compute_exposure_id(dayObs, seqNum, controller)
1241 else:
1242 raise ValueError("Unknown instrument {instrument}")
1245def computeCcdExposureId(instrument: str, exposureId: int, detector: int) -> int:
1246 instrument = instrument.lower()
1247 if instrument == "latiss":
1248 if detector != 0:
1249 raise ValueError("Invalid detector {detector} for LATISS")
1250 return lsst.obs.lsst.translators.LatissTranslator.compute_detector_exposure_id(exposureId, detector)
1251 elif instrument == "lsstcomcam":
1252 if detector < 0 or detector >= 9:
1253 raise ValueError("Invalid detector {detector} for LSSTComCam")
1254 return lsst.obs.lsst.translators.LsstComCamTranslator.compute_detector_exposure_id(
1255 exposureId, detector
1256 )
1257 elif instrument == "lsstcomcamsim":
1258 if detector < 0 or detector >= 9:
1259 raise ValueError("Invalid detector {detector} for LSSTComCamSim")
1260 return lsst.obs.lsst.translators.LsstComCamSimTranslator.compute_detector_exposure_id(
1261 exposureId, detector
1262 )
1263 elif instrument == "lsstcam":
1264 if detector < 0 or detector >= 205:
1265 raise ValueError("Invalid detector {detector} for LSSTCam")
1266 return lsst.obs.lsst.translators.LsstCamTranslator.compute_detector_exposure_id(exposureId, detector)
1267 else:
1268 raise ValueError("Unknown instrument {instrument}")
1271def getDetectorIds(instrumentName: str) -> list[int]:
1272 """Get a list of detector IDs for a given instrument.
1274 Parameters
1275 ----------
1276 instrumentName : `str`
1277 The name of the instrument.
1279 Returns
1280 -------
1281 detectorIds : `list` of `int`
1282 The list of detector IDs.
1283 """
1284 camera = getCameraFromInstrumentName(instrumentName)
1285 return [detector.getId() for detector in camera]
1288def getCameraFromInstrumentName(instrumentName: str) -> Camera:
1289 """Get the camera object given the instrument name (case insenstive).
1291 Parameters
1292 ----------
1293 instrumentName : `str`
1294 The name of the instrument, e.g. "LATISS" or "LSSTCam". Case
1295 insenstive.
1297 Returns
1298 -------
1299 camera: `lsst.afw.cameraGeom.Camera`
1300 The camera object corresponding to the instrument name.
1302 Raises
1303 ------
1304 ValueError
1305 If the instrument name is not supported.
1306 """
1308 _instrument = instrumentName.lower()
1310 match _instrument:
1311 case "lsstcam":
1312 camera = LsstCam.getCamera()
1313 case "lsstcomcam":
1314 camera = LsstComCam.getCamera()
1315 case "lsstcomcamsim":
1316 camera = LsstComCamSim.getCamera()
1317 case "latiss":
1318 camera = Latiss.getCamera()
1319 case _:
1320 raise ValueError(f"Unsupported instrument: {instrumentName}")
1321 return camera
1324def calcEclipticCoords(ra: float, dec: float) -> tuple[float, float]:
1325 """Get the ecliptic coordinates of the specified ra and dec.
1327 Transform J2000 (ra, dec), both in degrees, to
1328 IAU1976 Ecliptic coordinates (also returning degrees).
1330 Matches the results of:
1332 from astropy.coordinates import SkyCoord, HeliocentricEclipticIAU76
1333 from dataclasses import dataclass
1334 import astropy.units as u
1336 p = SkyCoord(ra=ra0*u.deg, dec=dec0*u.deg, distance=1*u.au, frame='hcrs')
1337 ecl = p.transform_to(HeliocentricEclipticIAU76)
1338 print(ecl.lon.value, ecl.lat.value)
1340 except that it's fast.
1342 Parameters
1343 ----------
1344 ra : `float`
1345 The right ascension, in degrees.
1346 dec : `float`
1347 The declination, in degrees.
1349 Returns
1350 -------
1351 lambda : `float`
1352 The ecliptic longitude in degrees.
1353 beta : `float`
1354 The ecliptic latitude in degrees.
1355 """
1357 ra, dec = np.deg2rad(ra), np.deg2rad(dec)
1359 # IAU 1976 obliquity
1360 epsilon = np.deg2rad(23.43929111111111)
1361 cos_eps, sin_eps = np.cos(epsilon), np.sin(epsilon)
1363 sra, cra = np.sin(ra), np.cos(ra)
1364 sdec, cdec = np.sin(dec), np.cos(dec)
1366 lambda_ = np.arctan2((sra * cos_eps + sdec / cdec * sin_eps), cra)
1367 beta = np.arcsin(sdec * cos_eps - cdec * sin_eps * sra)
1369 # normalize
1370 if lambda_ < 0:
1371 lambda_ += np.pi * 2
1373 return (np.rad2deg(lambda_), np.rad2deg(beta))
1376def getImageArray(
1377 inputData: np.ndarray | afwImage.Exposure | afwImage.Image | afwImage.MaskedImage,
1378) -> np.ndarray:
1379 """Get the image data from anything image-like.
1381 Parameters
1382 ----------
1383 inputData : `numpy.ndarray`, `lsst.afw.image.Exposure`,
1384 `lsst.afw.image.Image`, or `lsst.afw.image.MaskedImage`
1385 The input data.
1386 Returns
1387 -------
1388 imageData : `numpy.ndarray`
1389 The image data.
1390 Raises
1391 ------
1392 TypeError
1393 If the input data is not a numpy array, lsst.afw.image.Exposure,
1394 lsst.afw.image.Image, or lsst.afw.image.MaskedImage.
1395 """
1396 match inputData:
1397 case np.ndarray():
1398 imageData = inputData
1399 case afwImage.MaskedImage():
1400 imageData = inputData.image.array
1401 case afwImage.Image():
1402 imageData = inputData.array
1403 case afwImage.Exposure():
1404 imageData = inputData.image.array
1405 case _:
1406 raise TypeError(
1407 "This function accepts numpy array, lsst.afw.image.Exposure components."
1408 f" Got {type(inputData)}"
1409 )
1410 return imageData
1413def getSunAngle(time: Time | None = None) -> float:
1414 """Get the angle of the sun to the horizon at the specified time.
1416 If no time is specified, the current time is used. Positive numbers means
1417 the sun is above the horizon, negative means it is below.
1419 Returns
1420 -------
1421 sun_alt : `float`
1422 The angle of the sun to the horizon, in degrees.
1423 """
1424 if time is None:
1425 time = Time.now()
1426 sun_altaz = get_sun(time).transform_to(AltAz(obstime=time, location=SIMONYI_LOCATION))
1427 return sun_altaz.alt.deg
1430@dataclass(slots=True)
1431class RobustFitResult:
1432 slope: float
1433 intercept: float
1434 scatter: float
1435 slopePValue: float
1436 slopeStdErr: float
1437 slopeTValue: float
1438 interceptPValue: float
1439 interceptStdErr: float
1440 interceptTValue: float
1441 outlierMask: np.ndarray
1444class RobustFitter:
1445 """
1446 Robust linear fit using RANSAC + OLS, reusable across datasets.
1448 This class fits a robust linear trend to x and y using scikit-learn's
1449 RANSAC for inlier detection, followed by an OLS fit on inliers to compute
1450 slope, intercept, scatter, standard errors, t-values, and p-values.
1451 Configuration (e.g., minSamples) is set at init. Call `fit()` per dataset.
1452 Results are returned as a `RobustFitResult` and stored internally for later
1453 access (e.g., plotting).
1455 Parameters
1456 ----------
1457 minSamples : `float`, optional
1458 Minimum fraction of samples chosen randomly from the original data.
1459 """
1461 def __init__(self, *, minSamples: float = 0.2) -> None:
1462 self.minSamples = minSamples
1463 self._clearState()
1465 def _clearState(self) -> None:
1466 """Reset stored fit state."""
1467 self.model: Any = None
1468 self.olsModel: Any = None
1469 self.outlierMask: np.ndarray | None = None
1470 self.x: np.ndarray | None = None
1471 self.y: np.ndarray | None = None
1472 self.slope: float = np.nan
1473 self.intercept: float = np.nan
1474 self.ciSlope: tuple[float, float] = (np.nan, np.nan)
1475 self.ciIntercept: tuple[float, float] = (np.nan, np.nan)
1476 self.slopePValue = np.nan
1477 self.slopeStdErr = np.nan
1478 self.slopeTValue = np.nan
1479 self.scatter = np.nan
1481 @staticmethod
1482 def _defaultResidualThreshold(y: np.ndarray) -> float:
1483 """Compute default residual threshold from finite y values."""
1484 yFinite = np.asarray(y)[np.isfinite(y)]
1485 if yFinite.size == 0:
1486 raise ValueError("Cannot compute residual threshold: no finite y.")
1487 return 1.5 * float(mad_std(yFinite))
1489 def fit(
1490 self,
1491 x: np.ndarray,
1492 y: np.ndarray,
1493 *,
1494 residualThreshold: float | None = None,
1495 randomState: int = 42,
1496 ) -> RobustFitResult:
1497 """
1498 Fit robust line to (x, y) using RANSAC + OLS on inliers.
1500 Recomputes residual threshold each call if not provided.
1502 Parameters
1503 ----------
1504 x : array-like
1505 Independent variable values.
1506 y : array-like
1507 Dependent variable values.
1508 residualThreshold : float, optional
1509 Residual threshold for inlier detection. If None, computed
1510 as 1.5 × MAD of y.
1511 randomState : int, optional
1512 Random seed for RANSAC.
1514 Returns
1515 -------
1516 result : `RobustFitResult`
1517 Best-fit parameters and statistics from the OLS inlier fit.
1519 Raises
1520 ------
1521 ValueError
1522 If no finite x/y values are available for fitting.
1523 """
1524 self._clearState()
1526 x = np.asarray(x)
1527 y = np.asarray(y)
1529 finiteMask = np.isfinite(x) & np.isfinite(y)
1530 if not finiteMask.any():
1531 raise ValueError("No finite x/y values to fit.")
1533 X = x[finiteMask].reshape(-1, 1)
1534 yFit = y[finiteMask]
1536 if residualThreshold is None:
1537 residualThreshold = self._defaultResidualThreshold(yFit)
1539 ransac = RANSACRegressor(
1540 estimator=LinearRegression(),
1541 min_samples=self.minSamples,
1542 residual_threshold=residualThreshold,
1543 random_state=randomState,
1544 )
1545 ransac.fit(X, yFit)
1547 slope = float(ransac.estimator_.coef_[0])
1548 intercept = float(ransac.estimator_.intercept_)
1550 self.x = x
1551 self.y = y
1553 outlierMask = np.ones_like(x, dtype=bool)
1554 outlierMask[finiteMask] = ~ransac.inlier_mask_
1555 self.outlierMask = outlierMask
1557 xIn = x[~outlierMask]
1558 yIn = y[~outlierMask]
1559 XDesign = sm.add_constant(xIn)
1560 ols = sm.OLS(yIn, XDesign).fit()
1561 ci = ols.conf_int(alpha=0.32)
1562 ciIntercept = (float(ci[0][0]), float(ci[0][1]))
1563 ciSlope = (float(ci[1][0]), float(ci[1][1]))
1564 scatter = float(np.nanstd(ols.resid, ddof=1)) if ols.resid.size > 1 else np.nan
1566 self.model = ransac
1567 self.olsModel = ols
1568 self.slope = slope
1569 self.intercept = intercept
1570 self.ciSlope = ciSlope
1571 self.ciIntercept = ciIntercept
1572 self.slopePValue = float(ols.pvalues[1]) if np.isfinite(ols.pvalues[1]) else np.nan
1573 self.slopeStdErr = float(ols.bse[1]) if np.isfinite(ols.bse[1]) else np.nan
1574 self.slopeTValue = float(ols.tvalues[1]) if np.isfinite(ols.tvalues[1]) else np.nan
1575 self.scatter = scatter
1577 return RobustFitResult(
1578 slope=slope,
1579 intercept=intercept,
1580 scatter=scatter,
1581 slopePValue=self.slopePValue,
1582 slopeStdErr=self.slopeStdErr,
1583 slopeTValue=self.slopeTValue,
1584 interceptPValue=float(ols.pvalues[0]) if np.isfinite(ols.pvalues[0]) else np.nan,
1585 interceptStdErr=float(ols.bse[0]) if np.isfinite(ols.bse[0]) else np.nan,
1586 interceptTValue=float(ols.tvalues[0]) if np.isfinite(ols.tvalues[0]) else np.nan,
1587 outlierMask=self.outlierMask,
1588 )
1590 def reset(self) -> None:
1591 """Clear any stored fit state."""
1592 self._clearState()
1594 def plotBestFit(self, ax: matplotlib.axes.Axes, label=None, color=None, alphaBand=0.2, lw=2, nBins=5):
1595 """Plot the best fit line, confidence interval,
1596 and optionally scatter/binned data.
1598 Parameters
1599 ----------
1600 ax : matplotlib.axes.Axes, optional
1601 The axes to plot on. If None, uses current axes.
1602 label : str, optional
1603 Label for the best fit line.
1604 color : str or None, optional
1605 Color for the fit line and confidence band.
1606 alphaBand : float, optional
1607 Alpha transparency for the confidence interval band.
1608 lw : int, optional
1609 Line width for the best fit line.
1610 nBins : int, optional
1611 Number of bins for binned statistics.
1613 Returns
1614 -------
1615 ax : matplotlib.axes.Axes
1616 The axes with the plot.
1617 """
1618 # Handle case where self.x contains only NaNs
1619 if self.x is None or self.y is None or self.outlierMask is None or self.olsModel is None:
1620 raise RuntimeError("Fit must be called before plotting.")
1622 if np.isnan(self.x).all():
1623 raise ValueError("All x values are NaN; cannot plot best fit.")
1625 xx = np.linspace(np.nanmin(self.x), np.nanmax(self.x), 200)
1626 X_design = sm.add_constant(xx)
1628 pred = self.olsModel.get_prediction(X_design)
1629 summary_frame = pred.summary_frame(alpha=0.05) # 95% intervals
1631 mean = summary_frame["mean"]
1632 ci_lo = summary_frame["mean_ci_lower"]
1633 ci_hi = summary_frame["mean_ci_upper"]
1635 # Plot best fit
1636 ax.plot(xx, mean, color=color, label=label, lw=lw)
1637 # Plot confidence interval band
1638 ax.fill_between(xx, ci_lo, ci_hi, color=color, alpha=alphaBand)
1640 mask = self.outlierMask
1641 xin, yin = self.x[~mask], self.y[~mask] # inliers
1642 xout, _ = self.x[mask], self.y[mask] # outliers
1644 # Bin data in nbins
1645 if np.all(np.isnan(xout)):
1646 # Skip binning if all values are NaN
1647 pass
1648 else:
1649 bin_centers, means, stds, bin_counts = self.getBinnedData(xin, yin, nBins)
1651 # Plot binned means and stds
1652 ax.errorbar(
1653 bin_centers,
1654 means,
1655 color=color or "black",
1656 yerr=stds,
1657 fmt="o",
1658 capsize=4,
1659 markersize=6,
1660 alpha=0.7,
1661 zorder=10,
1662 )
1664 return ax
1666 @staticmethod
1667 def getBinnedData(
1668 x: np.ndarray, y: np.ndarray, nbins: int
1669 ) -> tuple[np.ndarray, list[float], list[float], list[int]]:
1670 """Get binned statistics for x and y."""
1671 bins = np.linspace(np.nanmin(x), np.nanmax(x), nbins + 1)
1672 bin_centers = 0.5 * (bins[:-1] + bins[1:])
1673 digitized = np.digitize(x, bins) - 1
1674 means = []
1675 stds = []
1676 bin_counts = []
1677 for i in range(nbins):
1678 bin_y = y[digitized == i]
1679 if len(bin_y) > 0:
1680 means.append(np.nanmedian(bin_y))
1681 if len(bin_y) > 1:
1682 stds.append(np.nanstd(bin_y, ddof=1) / np.sqrt(len(bin_y)))
1683 else:
1684 stds.append(np.nan)
1685 bin_counts.append(len(bin_y))
1686 else:
1687 means.append(np.nan)
1688 stds.append(np.nan)
1689 bin_counts.append(0)
1690 return bin_centers, means, stds, bin_counts