Coverage for python / lsst / summit / utils / astrometry / anet.py: 21%
184 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-25 09:03 +0000
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/>.
22import os
23import shutil
24import subprocess
25import tempfile
26import time
27import uuid
28import warnings
29from dataclasses import dataclass
30from functools import cached_property
31from typing import Any
33import numpy as np
34from astropy.io import fits
36import lsst.afw.geom as afwGeom
37import lsst.afw.image as afwImage
38import lsst.afw.table as afwTable
39import lsst.geom as geom
41from .utils import headerToWcs
43__all__ = ["AstrometryNetResult", "CommandLineSolver", "OnlineSolver"]
46@dataclass(frozen=True)
47class AstrometryNetResult:
48 """Minimal wrapper class to construct and return results from the command
49 line fitter.
51 Constructs a DM wcs from the output of the command line fitter, and
52 calculates the plate scale and astrometric scatter measurement in arcsec
53 and pixels.
55 Parameters
56 ----------
57 wcsFile : `str`
58 The path to the .wcs file from the fit.
59 corrFile : `str`, optional
60 The path to the .corr file from the fit.
61 """
63 wcsFile: str
64 corrFile: str | None = None
66 def __post_init__(self) -> None:
67 # touch these properties to ensure the files needed to calculate them
68 # are read immediately, in case they are deleted from temp
69 self.wcs
70 self.rmsErrorArsec
72 @cached_property
73 def wcs(self) -> afwGeom.SkyWcs:
74 with fits.open(self.wcsFile) as f:
75 header = f[0].header
76 return headerToWcs(header)
78 @cached_property
79 def plateScale(self) -> float:
80 return self.wcs.getPixelScale().asArcseconds()
82 @cached_property
83 def meanSqErr(self) -> float | None:
84 if not self.corrFile:
85 return None
87 try:
88 with fits.open(self.corrFile) as f:
89 data = f[1].data
91 meanSqErr = 0.0
92 count = 0
93 for i in range(data.shape[0]):
94 row = data[i]
95 count += 1
96 error = (row[0] - row[4]) ** 2 + (row[1] - row[5]) ** 2 # square error in pixels
97 error *= row[10] # multiply by weight
98 meanSqErr += error
99 meanSqErr /= count # divide by number of stars
100 return meanSqErr
101 except Exception as e:
102 print(f"Failed for calculate astrometric scatter: {repr(e)}")
103 return None
105 @cached_property
106 def rmsErrorPixels(self) -> float:
107 assert self.meanSqErr is not None, "No meanSqErr calculated, cannot calculate rmsErrorPixels."
108 return np.sqrt(self.meanSqErr)
110 @cached_property
111 def rmsErrorArsec(self) -> float:
112 return self.rmsErrorPixels * self.plateScale
115class CommandLineSolver:
116 """An interface for the solve-field command line tool from astrometry.net.
118 Parameters
119 ----------
120 indexFilePath : `str`
121 The path to the index files. Do not include the 4100 or 4200 etc. in
122 the path. This is selected automatically depending on the `isWideField`
123 flag when calling `run()`.
124 checkInParallel : `bool`, optional
125 Do the checks in parallel. Default is True.
126 timeout : `float`, optional
127 The timeout for the solve-field command. Default is 300 seconds.
128 binary : `str`, optional
129 The path to the solve-field binary. Default is 'solve-field', i.e. it
130 is assumed to be on the path.
131 """
133 def __init__(
134 self,
135 indexFilePath: str,
136 checkInParallel: bool = True,
137 timeout: float | int = 300,
138 binary: str = "solve-field",
139 fluxSlot: str = "base_CircularApertureFlux_3_0_instFlux",
140 ):
141 self.indexFilePath = indexFilePath
142 self.checkInParallel = checkInParallel
143 self.timeout = timeout
144 self.binary = binary
145 self.fluxSlot = fluxSlot
146 if not shutil.which(binary):
147 raise RuntimeError(
148 f"Could not find {binary} in path, please install 'solve-field' and either"
149 " put it on your PATH or specify the full path to it in the 'binary' argument"
150 )
152 def _writeConfigFile(self, wide: bool, useGaia: bool) -> str:
153 """Write a temporary config file for astrometry.net.
155 Parameters
156 ----------
157 wide : `bool`
158 Is this a wide field image? Used to select the 4100 vs 4200 dir in
159 the index file path. Ignored if ``useGaia`` is ``True``.
160 useGaia : `bool`
161 Use the 5200 Gaia catalog? If ``True``, ``wide`` is ignored.
163 Returns
164 -------
165 filename : `str`
166 The filename to which the config file was written.
167 """
168 fileSet = "4100" if wide else "4200"
169 fileSet = "5200/LITE" if useGaia else fileSet
170 indexFileDir = os.path.join(self.indexFilePath, fileSet)
171 if not os.path.isdir(indexFileDir):
172 raise RuntimeError(
173 f"No index files found at {self.indexFilePath}, in {indexFileDir} (you need a"
174 " 4100 dir for wide field and 4200 dir for narrow field images)."
175 )
177 lines = []
178 if self.checkInParallel:
179 lines.append("inparallel")
181 lines.append(f"cpulimit {self.timeout}")
182 lines.append(f"add_path {indexFileDir}")
183 lines.append("autoindex")
184 filename = tempfile.mktemp(suffix=".cfg")
185 with open(filename, "w") as f:
186 f.writelines(line + "\n" for line in lines)
187 return filename
189 def _writeFitsTable(self, sourceCat: afwTable.SourceCatalog) -> str:
190 """Write the source table to a FITS file and return the filename.
192 Parameters
193 ----------
194 sourceCat : `lsst.afw.table.SourceCatalog`
195 The source catalog to write to a FITS file for the solver.
197 Returns
198 -------
199 filename : `str`
200 The filename to which the catalog was written.
201 """
202 fluxArray = sourceCat[self.fluxSlot]
203 fluxFinite = np.logical_and(np.isfinite(fluxArray), fluxArray > 0)
204 fluxArray = fluxArray[fluxFinite]
205 indices = np.argsort(fluxArray)
206 x = sourceCat.getColumnView().getX()[fluxFinite]
207 y = sourceCat.getColumnView().getY()[fluxFinite]
208 fluxArray = fluxArray[indices][::-1] # brightest finite flux
209 xArray = x[indices][::-1]
210 yArray = y[indices][::-1]
211 x = fits.Column(name="X", format="D", array=xArray)
212 y = fits.Column(name="Y", format="D", array=yArray)
213 flux = fits.Column(name="FLUX", format="D", array=fluxArray)
214 print(f" of which {len(fluxArray)} made it into the fit")
215 hdu = fits.BinTableHDU.from_columns([flux, x, y])
217 filename = tempfile.mktemp(suffix=".fits")
218 hdu.writeto(filename)
219 return filename
221 # try to keep this call sig and the defaults as similar as possible
222 # to the run method on the OnlineSolver
223 def run(
224 self,
225 exp: afwImage.Exposure,
226 sourceCat: afwTable.SourceCatalog,
227 isWideField: bool,
228 *,
229 useGaia: bool = False,
230 percentageScaleError: float | int = 10,
231 radius: float | None = None,
232 silent: bool = True,
233 ) -> AstrometryNetResult | None:
234 """Get the astrometric solution for an image using astrometry.net using
235 the binary ``solve-field`` and a set of index files.
237 Parameters
238 ----------
239 exp : `lsst.afw.image.Exposure`
240 The input exposure. Only used for its wcs and its dimensions.
241 sourceCat : `lsst.afw.table.SourceCatalog`
242 The detected source catalog for the exposure. One produced by a
243 default run of CharacterizeImageTask is suitable.
244 isWideField : `bool`
245 Is this a wide field image? Used to select the correct index files.
246 Ignored if ``useGaia`` is ``True``.
247 useGaia : `bool`
248 Use the Gaia 5200/LITE index files? If set, ``isWideField`` is
249 ignored.
250 percentageScaleError : `float`, optional
251 The percentage scale error to allow in the astrometric solution.
252 radius : `float`, optional
253 The search radius from the nominal wcs in degrees.
254 silent : `bool`, optional
255 Swallow the output from the command line? The solver is *very*
256 chatty, so this is recommended.
258 Returns
259 -------
260 result : `AstrometryNetResult` or `None`
261 The result of the fit. If the fit was successful, the result will
262 contain a valid DM wcs, a scatter in arcseconds and a scatter in
263 pixels. If the fit failed, ``None`` is returned.
264 """
265 wcs = exp.getWcs()
266 if not wcs:
267 raise ValueError("No WCS in exposure")
269 configFile = self._writeConfigFile(wide=isWideField, useGaia=useGaia)
270 print(f"Fitting image with {len(sourceCat)} sources", end="")
271 fitsFile = self._writeFitsTable(sourceCat)
273 plateScale = wcs.getPixelScale().asArcseconds()
274 scaleMin = plateScale * (1 - percentageScaleError / 100)
275 scaleMax = plateScale * (1 + percentageScaleError / 100)
277 ra, dec = wcs.getSkyOrigin()
279 # do not use tempfile.TemporaryDirectory() because it must not exist,
280 # it is made by the solve-field binary and barfs if it exists already!
281 mainTempDir = tempfile.gettempdir()
282 tempDirSuffix = str(uuid.uuid1()).split("-")[0]
283 tempDir = os.path.join(mainTempDir, tempDirSuffix)
285 cmd = (
286 f"{self.binary} {fitsFile} " # the data
287 f"--width {exp.getWidth()} " # image dimensions
288 f"--height {exp.getHeight()} " # image dimensions
289 f"-3 {ra.asDegrees()} "
290 f"-4 {dec.asDegrees()} "
291 f"-5 {radius if radius else 180} "
292 "-X X -Y Y -v -z 2 -t 2 " # the parts of the bintable to use
293 f"--scale-low {scaleMin:.3f} " # the scale range
294 f"--scale-high {scaleMax:.3f} " # the scale range
295 f"--scale-units arcsecperpix "
296 f"--crpix-x {wcs.getPixelOrigin()[0]} " # set the pixel origin
297 f"--crpix-y {wcs.getPixelOrigin()[1]} " # set the pixel origin
298 f"--config {configFile} "
299 f"-D {tempDir} "
300 "--no-plots " # don't make plots
301 "--overwrite " # shouldn't matter as we're using temp files
302 )
304 t0 = time.time()
305 with open(os.devnull, "w") as devnull:
306 result = subprocess.run(cmd, shell=True, check=True, stdout=devnull if silent else None)
307 t1 = time.time()
309 if result.returncode == 0:
310 print(f"Fitting code ran in {(t1 - t0):.2f} seconds")
311 # output template is /tmpdirname/fitstempname + various suffixes
312 # for each obj
313 basename = os.path.basename(fitsFile).removesuffix(".fits")
314 outputTemplate = os.path.join(tempDir, basename)
315 wcsFile = outputTemplate + ".wcs"
316 corrFile = outputTemplate + ".corr"
318 if not os.path.exists(wcsFile):
319 print("but failed to find a solution.")
320 return None
322 fitResult = AstrometryNetResult(wcsFile, corrFile)
323 return fitResult
324 else:
325 print("Fit failed")
326 return None
329class OnlineSolver:
330 """A class to solve an image using the Astrometry.net online service."""
332 def __init__(self) -> None:
333 # import seems to spew warnings even if the required key is present
334 # so we swallow them, and raise on init if the key is missing
335 with warnings.catch_warnings():
336 warnings.simplefilter("ignore")
337 from astroquery.astrometry_net import AstrometryNet
339 self.apiKey = self.getApiKey() # raises if not present so do first
340 self.adn = AstrometryNet()
341 self.adn.api_key = self.apiKey
343 @staticmethod
344 def getApiKey() -> str:
345 """Get the astrometry.net API key if possible.
347 Raises a RuntimeError if it isn't found.
349 Returns
350 -------
351 apiKey : str
352 The astrometry.net API key, if present.
354 Raises
355 ------
356 RuntimeError
357 Raised if the ASTROMETRY_NET_API_KEY is not set.
358 """
359 try:
360 return os.environ["ASTROMETRY_NET_API_KEY"]
361 except KeyError as e:
362 msg = "No AstrometryNet API key found. Sign up and get one, set it to $ASTROMETRY_NET_API_KEY"
363 raise RuntimeError(msg) from e
365 # try to keep this call sig and the defaults as similar as possible
366 # to the run method on the CommandLineSolver
367 def run(
368 self,
369 exp: afwImage.Exposure,
370 sourceCat: afwTable.SourceCatalog,
371 *,
372 percentageScaleError: float | int = 10,
373 radius: float | None = None,
374 scaleEstimate: float | None = None,
375 ) -> dict[str, Any] | None:
376 """Get the astrometric solution for an image using the astrometry.net
377 online solver.
379 Parameters
380 ----------
381 exp : `lsst.afw.image.Exposure`
382 The input exposure. Only used for its wcs.
383 sourceCat : `lsst.afw.table.SourceCatalog`
384 The detected source catalog for the exposure. One produced by a
385 default run of CharacterizeImageTask is suitable.
386 percentageScaleError : `float`, optional
387 The percentage scale error to allow in the astrometric solution.
388 radius : `float`, optional
389 The search radius from the nominal wcs in degrees.
390 scaleEstimate : `float`, optional
391 An estimate of the scale in arcseconds per pixel. Only used if
392 (and required when) the exposure has no wcs.
394 Returns
395 -------
396 result : `dict` or `None`
397 The results of the fit, with the following keys, or ``None`` if
398 the fit failed:
399 ``nominalRa`` : `lsst.geom.Angle`
400 The nominal ra from the exposure's boresight.
401 ``nominalDec`` : `lsst.geom.Angle`
402 The nominal dec from the exposure's boresight.
403 ``calculatedRa`` : `lsst.geom.Angle`
404 The fitted ra.
405 ``calculatedDec`` : `lsst.geom.Angle`
406 The fitted dec.
407 ``deltaRa`` : `lsst.geom.Angle`,
408 The change in ra, as an Angle.
409 ``deltaDec`` : `lsst.geom.Angle`,
410 The change in dec, as an Angle.
411 ``deltaRaArcsec`` : `float``
412 The change in ra in arcseconds, as a float.
413 ``deltaDecArcsec`` : `float`
414 The change in dec in arcseconds, as a float.
415 ``astrometry_net_wcs_header`` : `dict`
416 The fitted wcs, as a header dict.
417 """
418 nominalWcs = exp.getWcs()
419 if nominalWcs is not None:
420 ra, dec = nominalWcs.getSkyOrigin()
421 scaleEstimate = nominalWcs.getPixelScale().asArcseconds()
422 else:
423 print("Trying to process image with None wcs - good luck!")
424 vi = exp.getInfo().getVisitInfo()
425 ra, dec = vi.boresightRaDec
426 if np.isnan(ra.asDegrees()) or np.isnan(dec.asDegrees()):
427 raise RuntimeError("Exposure has no wcs and did not find nominal ra/dec in visitInfo")
429 if not scaleEstimate: # must either have a wcs or provide via kwarg
430 raise RuntimeError("Got no kwarg for scaleEstimate and failed to find one in the nominal wcs.")
432 image_height, image_width = exp.image.array.shape
433 scale_units = "arcsecperpix"
434 scale_type = "ev" # ev means submit estimate and % error
435 scale_err = percentageScaleError # error as percentage
436 center_ra = ra.asDegrees()
437 center_dec = dec.asDegrees()
438 radius = radius if radius else 180 # degrees
439 try:
440 wcs_header = self.adn.solve_from_source_list(
441 sourceCat["base_SdssCentroid_x"],
442 sourceCat["base_SdssCentroid_y"],
443 image_width,
444 image_height,
445 scale_units=scale_units,
446 scale_type=scale_type,
447 scale_est=scaleEstimate,
448 scale_err=scale_err,
449 center_ra=center_ra,
450 center_dec=center_dec,
451 radius=radius,
452 crpix_center=True, # the CRPIX is always the center
453 solve_timeout=240,
454 )
455 except RuntimeError:
456 print("Failed to find a solution")
457 return None
459 print("Finished solving!")
461 nominalRa, nominalDec = exp.getInfo().getVisitInfo().getBoresightRaDec()
463 if "CRVAL1" not in wcs_header:
464 raise RuntimeError("Astrometric fit failed.")
465 calculatedRa = geom.Angle(wcs_header["CRVAL1"], geom.degrees)
466 calculatedDec = geom.Angle(wcs_header["CRVAL2"], geom.degrees)
468 deltaRa = geom.Angle(wcs_header["CRVAL1"] - nominalRa.asDegrees(), geom.degrees)
469 deltaDec = geom.Angle(wcs_header["CRVAL2"] - nominalDec.asDegrees(), geom.degrees)
471 # TODO: DM-37213 change this to return an AstrometryNetResult class
472 # like the CommandLineSolver does.
474 result = {
475 "nominalRa": nominalRa,
476 "nominalDec": nominalDec,
477 "calculatedRa": calculatedRa,
478 "calculatedDec": calculatedDec,
479 "deltaRa": deltaRa,
480 "deltaDec": deltaDec,
481 "deltaRaArcsec": deltaRa.asArcseconds(),
482 "deltaDecArcsec": deltaDec.asArcseconds(),
483 "astrometry_net_wcs_header": wcs_header,
484 "nSources": len(sourceCat),
485 }
487 return result