Coverage for python/lsst/analysis/tools/actions/keyedData/calcDistances.py: 80%
120 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-11 09:02 +0000
1# This file is part of analysis_tools.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21__all__ = ("CalcRelativeDistances",)
23import astropy.units as u
24import esutil
25import numpy as np
26from smatch import Matcher
28from lsst.pex.config import Field
30from ...interfaces import KeyedData, KeyedDataAction, KeyedDataSchema, Vector
33class CalcRelativeDistances(KeyedDataAction):
34 """Calculate relative distances in a matched catalog.
36 Given a catalog of matched sources from multiple visits, this finds all
37 pairs of objects at a given separation, then calculates the separation of
38 their component source measurements from the individual visits. The RMS of
39 these is used to calculate the astrometric relative repeatability metric,
40 AMx, while the overall distribution of separations is used to compute the
41 ADx and AFx metrics.
42 """
44 groupKey = Field[str](doc="Column key to use for forming groups", default="obj_index")
45 visitKey = Field[str](doc="Column key to use for matching visits", default="visit")
46 raKey = Field[str](doc="RA column key", default="coord_ra")
47 decKey = Field[str](doc="Dec column key", default="coord_dec")
48 annulus = Field[float](doc="Radial distance of the annulus in arcmin", default=5.0)
49 width = Field[float](doc="Width of annulus in arcmin", default=2.0)
50 threshAD = Field[float](doc="Threshold in mas for AFx calculation.", default=20.0)
51 threshAF = Field[float](
52 doc="Percentile of differences that can vary by more than threshAD.", default=10.0
53 )
54 maxPairs = Field[int](
55 doc="Maximum number of pairs to use; downsample otherwise.",
56 default=100_000,
57 )
58 bruteForceRadiusThreshold = Field[float](
59 doc="Matching radius at which to switch to brute-force matching (deg).",
60 default=0.1,
61 )
62 bruteForceMaxMatch = Field[int](
63 doc="Number of objects to sample for brute-force matching.",
64 default=10_000,
65 )
66 randomSeed = Field[int](
67 doc="Random seed to use when downsampling.",
68 default=12345,
69 )
71 def getInputSchema(self) -> KeyedDataSchema:
72 return (
73 (self.groupKey, Vector),
74 (self.raKey, Vector),
75 (self.decKey, Vector),
76 (self.visitKey, Vector),
77 )
79 def __call__(self, data: KeyedData, **kwargs) -> KeyedData:
80 """Run the calculation.
82 Parameters
83 ----------
84 data : KeyedData
85 Catalog of data including coordinate, visit, and object group
86 information.
88 Returns
89 -------
90 distanceParams : `dict`
91 Dictionary of the calculated arrays and metrics with the following
92 keys:
94 - ``rmsDistances`` : Per-object rms of separations (`np.array`).
95 - ``separationResiduals`` : All separations minus per-object median
96 (`np.array`)
97 - ``AMx`` : AMx metric (`float`).
98 - ``ADx`` : ADx metric (`float`).
99 - ``AFx`` : AFx metric (`float`).
100 """
101 distanceParams = {
102 "rmsDistances": np.array([]),
103 "separationResiduals": np.array([]),
104 "AMx": np.nan,
105 "ADx": np.nan,
106 "AFx": np.nan,
107 "nPairs": 0,
108 }
110 if len(data[self.groupKey]) == 0: 110 ↛ 111line 110 didn't jump to line 111 because the condition on line 110 was never true
111 return distanceParams
113 rng = np.random.RandomState(seed=self.randomSeed)
115 _, groupId = np.unique(data[self.groupKey], return_inverse=True)
116 nObj = groupId.max() + 1
118 # Compute the meanRa/meanDec.
119 meanRa = np.zeros(nObj)
120 meanDec = np.zeros_like(meanRa)
121 nObs = np.zeros_like(meanRa, dtype=np.int64)
123 # Check if tract is overlapping ra=0 and rotate if so.
124 # We assume a tract is smaller than 60x60 degrees.
125 rotation = 0.0
126 if np.max(data[self.raKey]) > 330.0 and np.min(data[self.raKey]) < 30.0: 126 ↛ 127line 126 didn't jump to line 127 because the condition on line 126 was never true
127 rotation = 180.0
128 raRotated = np.array(data[self.raKey]) - rotation
129 else:
130 raRotated = np.array(data[self.raKey])
132 np.add.at(meanRa, groupId, raRotated)
133 np.add.at(meanDec, groupId, np.array(data[self.decKey]))
134 np.add.at(nObs, groupId, 1)
136 meanRa /= nObs
137 meanDec /= nObs
138 meanRa += rotation
140 D = (self.annulus * u.arcmin).to_value(u.degree)
141 width = (self.width * u.arcmin).to_value(u.degree)
142 annulus = D + (width / 2) * np.array([-1, +1])
144 if annulus[1] >= self.bruteForceRadiusThreshold: 144 ↛ 146line 144 didn't jump to line 146 because the condition on line 144 was never true
145 # Use brute-force matching.
146 if len(meanRa) > self.bruteForceMaxMatch:
147 bruteForceSelection = rng.choice(len(meanRa), size=self.bruteForceMaxMatch, replace=False)
148 else:
149 bruteForceSelection = np.arange(len(meanRa))
151 i1 = np.tile(bruteForceSelection, len(bruteForceSelection))
152 i2 = np.repeat(bruteForceSelection, len(bruteForceSelection))
154 d = np.rad2deg(
155 sphDist(
156 np.deg2rad(meanRa)[i1],
157 np.deg2rad(meanDec)[i1],
158 np.deg2rad(meanRa)[i2],
159 np.deg2rad(meanDec)[i2],
160 ),
161 )
162 else:
163 with Matcher(meanRa, meanDec) as m:
164 idx, i1, i2, d = m.query_self(annulus[1], return_indices=True)
166 inAnnulus = (d > annulus[0]) & (d < annulus[1])
167 i1 = i1[inAnnulus]
168 i2 = i2[inAnnulus]
170 if len(i1) == 0: 170 ↛ 171line 170 didn't jump to line 171 because the condition on line 170 was never true
171 return distanceParams
173 if len(i1) > self.maxPairs: 173 ↛ 175line 173 didn't jump to line 175 because the condition on line 173 was never true
174 # Downsample the pairs.
175 selection = rng.choice(len(i1), size=self.maxPairs, replace=False)
176 i1 = i1[selection]
177 i2 = i2[selection]
179 # Match groups and get indices.
180 h, rev = esutil.stat.histogram(groupId, rev=True)
182 # Match together pairs that have the same visit.
183 # It is unfortunate that this requires a loop, but it is not slow.
184 # After this matching we have a set of matchedObsInd1/matchedObsInd2
185 # that are all individual observations that are in the annulus and
186 # share a visit. The matchedPairInd groups all the paired observations
187 # of a given pair.
188 matchedObsInd1 = []
189 matchedObsInd2 = []
190 matchedPairInd = []
191 for ind in range(len(i1)):
192 objInd1 = i1[ind]
193 objInd2 = i2[ind]
194 obsInd1 = rev[rev[objInd1] : rev[objInd1 + 1]]
195 obsInd2 = rev[rev[objInd2] : rev[objInd2 + 1]]
196 a, b = esutil.numpy_util.match(data[self.visitKey][obsInd1], data[self.visitKey][obsInd2])
197 matchedObsInd1.append(obsInd1[a])
198 matchedObsInd2.append(obsInd2[b])
199 matchedPairInd.append(np.full(len(a), ind))
201 matchedObsInd1 = np.concatenate(matchedObsInd1)
202 matchedObsInd2 = np.concatenate(matchedObsInd2)
203 matchedPairInd = np.concatenate(matchedPairInd)
205 separations = sphDist(
206 np.deg2rad(np.array(data[self.raKey][matchedObsInd1])),
207 np.deg2rad(np.array(data[self.decKey][matchedObsInd1])),
208 np.deg2rad(np.array(data[self.raKey][matchedObsInd2])),
209 np.deg2rad(np.array(data[self.decKey][matchedObsInd2])),
210 )
212 # Compute the mean from the ragged array of pairs by
213 # using np.add.at to sum numerator and denominator.
214 sepMean = np.zeros(len(i1))
215 nSep = np.zeros_like(sepMean, dtype=np.int32)
216 np.add.at(sepMean, matchedPairInd, separations)
217 np.add.at(nSep, matchedPairInd, 1)
218 good = nSep > 1
219 sepMean[good] /= nSep[good]
220 sepMean[~good] = np.nan
222 # There are no good pairs, so return the default.
223 if good.sum() == 0: 223 ↛ 224line 223 didn't jump to line 224 because the condition on line 223 was never true
224 return distanceParams
226 # Compute the stdev with sqrt(sum((sep - mean(sep))**2.)/(nsep - 1))
227 sepStd = np.zeros_like(sepMean)
228 np.add.at(
229 sepStd,
230 matchedPairInd,
231 (separations - sepMean[matchedPairInd]) ** 2.0,
232 )
233 sepStd[good] = np.sqrt(sepStd[good] / (nSep[good] - 1))
234 rmsDistances = sepStd[good]
236 # Need sepResiduals, but only when nSep is > 2.
237 bad2 = nSep <= 2
238 sepMean[bad2] = np.nan
239 sepResiduals = separations - sepMean[matchedPairInd]
240 sepResiduals = sepResiduals[np.isfinite(sepResiduals)]
242 # This is always going to be valid because we checked the number
243 # of good pairs above.
244 AMx = (np.median(rmsDistances) * u.radian).to(u.marcsec)
246 # Because there is a more stringent selection for sepResiduals,
247 # we need to check that we have enough to compute the metrics.
248 if len(sepResiduals) <= 1: 248 ↛ 249line 248 didn't jump to line 249 because the condition on line 248 was never true
249 AFx = np.nan * u.percent
250 ADx = np.nan * u.marcsec
251 absDiffSeparations = np.array([]) * u.marcsec
252 else:
253 absDiffSeparations = (abs(sepResiduals - np.median(sepResiduals)) * u.radian).to(u.marcsec)
254 afThreshhold = 100.0 - self.threshAF
255 ADx = np.percentile(absDiffSeparations, afThreshhold)
256 AFx = 100 * np.mean(np.abs(absDiffSeparations) > self.threshAD * u.marcsec) * u.percent
258 distanceParams["rmsDistances"] = (rmsDistances * u.radian).to(u.marcsec).value
259 distanceParams["separationResiduals"] = absDiffSeparations.value
260 distanceParams["AMx"] = AMx.value
261 distanceParams["ADx"] = ADx.value
262 distanceParams["AFx"] = AFx.value
263 distanceParams["nPairs"] = len(rmsDistances)
265 return distanceParams
268def sphDist(ra_mean, dec_mean, ra, dec):
269 """Calculate distance on the surface of a unit sphere.
271 This function was borrowed from faro.
273 Parameters
274 ----------
275 ra_mean : `float`
276 Mean RA in radians.
277 dec_mean : `float`
278 Mean Dec in radians.
279 ra : `numpy.array` [`float`]
280 Array of RA in radians.
281 dec : `numpy.array` [`float`]
282 Array of Dec in radians.
284 Notes
285 -----
286 Uses the Haversine formula to preserve accuracy at small angles.
287 Law of cosines approach doesn't work well for the typically very small
288 differences that we're looking at here.
289 """
290 # Haversine
291 dra = ra - ra_mean
292 ddec = dec - dec_mean
293 a = np.square(np.sin(ddec / 2)) + np.cos(dec_mean) * np.cos(dec) * np.square(np.sin(dra / 2))
294 dist = 2 * np.arcsin(np.sqrt(a))
296 # This is what the law of cosines would look like
297 # dist = np.arccos(np.sin(dec1)*np.sin(dec2) +
298 # np.cos(dec1)*np.cos(dec2)*np.cos(ra1 - ra2))
300 # This will also work, but must run separately for each element
301 # whereas the numpy version will run on either scalars or arrays:
302 # sp1 = geom.SpherePoint(ra1, dec1, geom.radians)
303 # sp2 = geom.SpherePoint(ra2, dec2, geom.radians)
304 # return sp1.separation(sp2).asRadians()
306 return dist