"""
Astronomical pointing correction and plate solving utilities.
This module provides comprehensive tools for correcting telescope pointing errors
by analyzing astronomical images and comparing detected stars with catalog data.
The primary workflow involves star detection, world coordinate system (WCS)
computation, and pointing offset calculation.
Key Features:
- Multi-scale star detection algorithms for robust identification in noisy images
- Gaia star catalog integration for precise astrometric reference
- World Coordinate System (WCS) computation using the twirl library
- Automatic pointing correction calculation and validation
- Support for FITS file processing with metadata extraction
- Background subtraction and image cleaning utilities
The module supports two main use cases:
1. Real-time pointing correction during observations
2. Post-processing analysis of astronomical images
Typical Workflow:
1. Load or generate an astronomical image
2. Detect stars using multi-scale algorithms
3. Query Gaia catalog for reference stars in the field
4. Compute WCS transformation between image and sky coordinates
5. Calculate pointing correction based on target vs actual position
6. Apply corrections to telescope pointing
Example:
.. code-block:: yaml
# Process a FITS file for pointing correction
pointing_correction, image_star_mapping, num_stars = (
calculate_pointing_correction_from_fits(
"observation.fits",
target_ra=150.25,
target_dec=2.18
)
)
# Get the pointing offset
ra_offset = pointing_correction.offset_ra
dec_offset = pointing_correction.offset_dec
print(f"Pointing error: RA={ra_offset:.3f}°, Dec={dec_offset:.3f}°")
"""
import logging
import sqlite3
from dataclasses import dataclass
from datetime import datetime
from pathlib import Path
from typing import Optional, Tuple, Union
import astropy.units as u
import cabaret
import numpy as np
import pandas as pd
import twirl
from astropy.coordinates import SkyCoord
from astropy.io import fits
from astropy.stats import sigma_clipped_stats
from astropy.units import Quantity
from astropy.wcs.utils import WCS, pixel_to_skycoord, proj_plane_pixel_scales
from matplotlib import pyplot as plt
from photutils.detection import DAOStarFinder
from astra.config import Config
from astra.utils import optics
from astra.utils.image import clean_image
logger = logging.getLogger(__name__)
[docs]
def calculate_pointing_correction_from_fits(
filepath: str | Path | None,
dark_frame: str | Path | np.ndarray | None = None,
target_ra: float | None = None,
target_dec: float | None = None,
filter_band: Optional[str] = None,
use_local_db: bool = False,
logger: Optional[logging.Logger] = None,
):
"""
Calculate pointing correction from a FITS file.
Performs the complete plate solving workflow from a FITS file, including
metadata extraction, image cleaning, star detection, and pointing correction.
Parameters:
filepath (str or Path): Path to the FITS file.
dark_frame (str, Path or np.ndarray, optional): Dark frame for
subtraction, as a path or as image data. Defaults to None.
target_ra (float, optional): Target right ascension in degrees. If not
provided, it's read from FITS header. Defaults to None.
target_dec (float, optional): Target declination in degrees. If not
provided, it's read from FITS header. Defaults to None.
filter_band (str, optional): Filter band used for observation. Defaults to None.
use_local_db (bool, optional): Whether to use the local Gaia database. Defaults to False.
logger (logging.Logger, optional): Logger for progress messages, for
example the observatory logger so the messages reach the UI.
Defaults to the module logger.
Returns:
Tuple[PointingCorrection, ImageStarMapping, int]: A tuple containing:
- PointingCorrection: Object with computed pointing correction
- ImageStarMapping: Mapping of image stars to Gaia stars
- int: Number of stars used from the image for plate solving
Raises:
Exception: If insufficient stars are detected for plate solving,
or if the offset is larger than the field of view.
PlateSolveError: If no solution passes the checks in
``_verify_plate_solve``.
"""
image, header = _read_fits_file(filepath)
if dark_frame is not None:
if not isinstance(dark_frame, np.ndarray):
dark_frame, _ = _read_fits_file(dark_frame)
image = image - dark_frame
if target_ra is None:
target_ra = header["RA"]
if target_dec is None:
target_dec = header["DEC"]
plate_scale, dateobs = _extract_plate_scale_and_dateobs(header)
return calculate_pointing_correction_from_image(
image,
target_ra=target_ra,
target_dec=target_dec,
dateobs=dateobs,
plate_scale=plate_scale,
filter_band=filter_band,
use_local_db=use_local_db,
logger=logger,
)
[docs]
def calculate_pointing_correction_from_image(
image: np.ndarray,
target_ra: float,
target_dec: float,
dateobs: datetime,
plate_scale: float,
filter_band: Optional[str] = None,
use_local_db: bool = False,
logger: Optional[logging.Logger] = None,
):
"""
Calculate pointing correction from an image array.
Performs the complete plate solving workflow: image cleaning, star detection,
Gaia catalog matching, WCS computation, and pointing correction calculation.
Parameters:
image (np.ndarray): The 2D astronomical image data.
target_ra (float): Target right ascension in degrees.
target_dec (float): Target declination in degrees.
dateobs (datetime): Observation date for proper motion corrections.
plate_scale (float): Image plate scale in degrees per pixel.
filter_band (str, optional): Filter band used for observation.
use_local_db (bool, optional): Whether to use the local Gaia database. Defaults to False.
logger (logging.Logger, optional): Logger for progress messages, for
example the observatory logger so the messages reach the UI.
Defaults to the module logger.
Returns:
Tuple[PointingCorrection, ImageStarMapping, int]: A tuple containing:
- PointingCorrection: Object with computed pointing correction
- ImageStarMapping: Mapping of image stars to Gaia stars
- int: Number of stars used from the image for plate solving
Raises:
Exception: If insufficient stars are detected for plate solving,
or if the offset is larger than the field of view.
PlateSolveError: If no solution passes the checks in
``_verify_plate_solve``.
"""
logger = logging.getLogger(__name__) if logger is None else logger
image_clean = clean_image(image)
# assume 2" FWHM
stars_in_image = find_stars(
image_clean,
threshold=7,
fwhm=optics.fwhm_pixels(
plate_scale_deg_per_pixel=plate_scale, seeing_arcsec=2.0
),
)
# Limit number of stars and gaia stars to use for plate solve
stars_in_image_used = min(len(stars_in_image), 20)
if stars_in_image_used < 4:
raise Exception("Not enough stars detected to plate solve")
stars_in_image = stars_in_image[0:stars_in_image_used]
logger.debug(f"Detected {len(stars_in_image)} stars in image for plate solving")
gaia_star_coordinates = _get_gaia_star_coordinates(
target_ra,
target_dec,
image_clean,
dateobs,
plate_scale,
filter_band=filter_band,
fov_scale=1.25,
limit=int(2 * stars_in_image_used),
use_local_db=use_local_db,
logger=logger,
)
logger.debug(f"Queried {len(gaia_star_coordinates)} Gaia stars for plate solving")
image_star_mapping = ImageStarMapping.from_gaia_coordinates(
stars_in_image, gaia_star_coordinates, plate_scale=plate_scale, logger=logger
)
plating_ra, plating_dec = image_star_mapping.get_plating_center(
image_shape=image_clean.shape
)
pointing_correction = PointingCorrection(
target_ra=target_ra,
target_dec=target_dec,
plating_ra=plating_ra,
plating_dec=plating_dec,
)
_verify_offset_within_fov(pointing_correction, plate_scale, image_clean.shape)
return pointing_correction, image_star_mapping, stars_in_image_used
[docs]
def find_stars(
data: np.ndarray,
threshold: float = 5.0,
fwhm: float = 3.0,
) -> np.ndarray:
"""
Find stars using DAOStarFinder algorithm.
Uses the photutils DAOStarFinder algorithm to detect point sources in
astronomical images. The function performs background subtraction and
returns star coordinates sorted by brightness.
Parameters:
data (np.ndarray): The 2D image data array.
threshold (float, optional): Detection threshold in units of background
standard deviation. Higher values detect fewer, brighter stars.
Defaults to 5.0.
fwhm (float, optional): Expected Full Width at Half Maximum of stars
in pixels. Should match the typical seeing conditions. Defaults to 3.0.
Returns:
np.ndarray: Array of detected star coordinates sorted by brightness.
Shape is (N, 2) where N is the number of stars, and each row is (x, y).
Returns an empty array if no stars are found.
"""
# Calculate background statistics
mean, median, std = sigma_clipped_stats(data, sigma=3.0)
# Use DAOStarFinder for star detection
dao_find = DAOStarFinder(
fwhm=fwhm,
threshold=threshold * std,
exclude_border=True,
min_separation=2 * fwhm,
)
dao_sources = dao_find(data)
if dao_sources is None or len(dao_sources) == 0:
return np.array([]).reshape(0, 2)
# Sort sources by flux (brightness) in descending order
sorted_indices = np.argsort(dao_sources["flux"])[::-1]
dao_sources = dao_sources[sorted_indices]
# Filter sources based on peak value
dao_sources = dao_sources[dao_sources["peak"] > mean + threshold * std]
# Convert to (x, y) coordinates
coordinates = np.column_stack(
[dao_sources["x_centroid"], dao_sources["y_centroid"]]
)
return coordinates
[docs]
class PlateSolveError(Exception):
"""Raised when twirl finds no solution, or no solution passes the checks."""
[docs]
@dataclass
class PointingCorrection:
"""
Class to store the pointing correction between the desired target center and the plating center.
Stores target and actual telescope pointing coordinates and provides methods
to calculate offsets, angular separations, and proxy coordinates for correction.
The proxy coordinates represent where the telescope should point to achieve
the desired target position.
Properties:
offset_ra: RA offset (plating - target) in degrees :no-index:
offset_dec: Dec offset (plating - target) in degrees :no-index:
angular_separation: Angular distance between target and actual position :no-index:
proxy_ra: RA coordinate to point telescope to reach target :no-index:
proxy_dec: Dec coordinate to point telescope to reach target :no-index:
Example:
>>> from astra.pointer import PointingCorrection
>>> pointing_correction = PointingCorrection(
... target_ra=10.685, target_dec=41.269,
... plating_ra=10.68471, plating_dec=41.26917
... )
>>> print(f"RA offset: {pointing_correction.offset_ra:.4f} degrees")
"""
"""The right ascension of the target center in degrees."""
target_ra: float
"""The declination of the target center in degrees."""
target_dec: float
"""The right ascension of the plating center in degrees."""
plating_ra: float
"""The declination of the plating center in degrees."""
plating_dec: float
@property
def offset_ra(self):
"""
Calculate the RA offset between plating and target coordinates.
Returns:
float: RA offset (plating_ra - target_ra) in degrees. Positive values
indicate the telescope pointed east of the target.
"""
return self.plating_ra - self.target_ra
@property
def offset_dec(self):
"""
Calculate the Dec offset between plating and target coordinates.
Returns:
float: Dec offset (plating_dec - target_dec) in degrees. Positive values
indicate the telescope pointed north of the target.
"""
return self.plating_dec - self.target_dec
@property
def angular_separation(self) -> float:
"""
Calculate the angular separation between target and plating coordinates.
Returns:
float: Angular distance between target and actual position in degrees.
Always positive, represents the magnitude of the pointing error.
"""
target_coord = SkyCoord(ra=self.target_ra, dec=self.target_dec, unit="deg")
plating_coord = SkyCoord(ra=self.plating_ra, dec=self.plating_dec, unit="deg")
return target_coord.separation(plating_coord).deg
@property
def proxy_ra(self):
"""
Calculate the RA coordinate to point the telescope to reach the target RA.
This assumes that the offset is independent of the telescope's position in the sky.
The proxy coordinate compensates for the systematic pointing error.
Returns:
float: RA coordinate in degrees where the telescope should point to
effectively arrive at the target RA.
See Also:
ConstantDistorter: Verify sign convention of correction.
"""
return self.target_ra - self.offset_ra
@property
def proxy_dec(self):
"""
Calculate the Dec coordinate to point the telescope to reach the target Dec.
This assumes that the offset is independent of the telescope's position in the sky.
The proxy coordinate compensates for the systematic pointing error.
Returns:
float: Dec coordinate in degrees where the telescope should point to
effectively arrive at the target Dec.
See Also:
ConstantDistorter: Verify sign convention of correction.
"""
return self.target_dec - self.offset_dec
def __repr__(self):
return (
"PointingCorrection("
f"target_ra={self.target_ra}, target_dec={self.target_dec}, "
f" plating_ra={self.plating_ra}, plating_dec={self.plating_dec})"
)
[docs]
@dataclass
class ImageStarMapping:
"""
A class to handle the mapping of stars detected in an image to their corresponding
Gaia star coordinates using World Coordinate System (WCS) transformations.
Manages the relationship between pixel coordinates of detected stars and their
celestial coordinates from the Gaia catalog. Provides methods for coordinate
transformations, star matching, and validation of the plate solving process.
"""
"""The World Coordinate System object used for mapping celestial coordinates to pixel coordinates."""
wcs: WCS
"""An array of detected star coordinates in the image, represented in pixel format (x, y)."""
stars_in_image: np.ndarray
"""An array of Gaia star coordinates projected into the image's pixel space using the WCS transformation."""
gaia_stars_in_image: np.ndarray
[docs]
@classmethod
def from_gaia_coordinates(
cls,
stars_in_image: np.ndarray,
gaia_stars: np.ndarray,
plate_scale: Optional[float] = None,
logger: Optional[logging.Logger] = None,
):
"""
Create an ImageStarMapping instance from detected stars and Gaia coordinates.
Computes the World Coordinate System (WCS) transformation using the twirl
library to map between pixel and sky coordinates.
If ``plate_scale`` is given, twirl first matches triangles of stars,
which is fast. If that gives no solution, raises, or the solution
fails ``_verify_plate_solve``, it tries once more with quads of stars
(the twirl defaults). Triangles can give wrong solutions that only the
checks catch, so without ``plate_scale`` only quads are used and the
solution is not checked.
Parameters:
stars_in_image (np.ndarray): Array of detected star coordinates in pixels (x, y).
gaia_stars (np.ndarray): Array of Gaia star coordinates in degrees (ra, dec).
plate_scale (float, optional): Expected plate scale in degrees per
pixel. If given, each solution must pass ``_verify_plate_solve``.
logger (logging.Logger, optional): Logger for the result of each
attempt. Defaults to the module logger.
Returns:
ImageStarMapping: New instance with computed WCS and transformed coordinates.
Raises:
PlateSolveError: If no attempt gives a solution that passes the
checks.
"""
logger = logging.getLogger(__name__) if logger is None else logger
if plate_scale is None:
attempts = {"quads": {}}
else:
attempts = {
"triangles": {"asterism": 3, "min_match": 0.5},
"quads": {},
}
failures = []
for name, kwargs in attempts.items():
try:
wcs = twirl.compute_wcs(stars_in_image, gaia_stars, **kwargs)
except Exception as e:
# e.g. fit_wcs_from_points raises on a bad match
failures.append(f"{name}: {e}")
logger.info(f"Plate solve with {name} raised: {e}")
continue
if wcs is None:
failures.append(
f"{name}: no match between the {len(stars_in_image)} "
f"image stars and the {len(gaia_stars)} Gaia stars"
)
logger.info(f"Plate solve with {name} found no match")
continue
gaia_stars_in_image = np.array(
SkyCoord(gaia_stars, unit="deg").to_pixel(wcs)
).T
mapping = cls(wcs, stars_in_image, gaia_stars_in_image)
if plate_scale is None:
return mapping
try:
summary = _verify_plate_solve(mapping, plate_scale)
except PlateSolveError as e:
failures.append(f"{name}: {e}")
logger.info(f"Plate solve with {name} rejected: {e}")
continue
fallback = "quads fallback" if failures else "no fallback"
logger.info(f"Plate solved with {name} ({fallback}): {summary}")
return mapping
raise PlateSolveError("Plate solve failed. " + "; ".join(failures))
[docs]
def get_plating_center(self, image_shape: Tuple[int, int]) -> Tuple[float, float]:
"""
Calculate the sky coordinates of the image center (plating center).
Parameters:
image_shape (Tuple[int, int]): Shape of the image as (height, width).
Returns:
Tuple[float, float]: RA and Dec coordinates of the image center in degrees.
"""
plating_center = pixel_to_skycoord(
image_shape[1] / 2, image_shape[0] / 2, self.wcs
)
return plating_center.ra.deg, plating_center.dec.deg
[docs]
def skycoord_to_pixels(self, ra: float, dec: float) -> Tuple[float, float]:
"""
Convert sky coordinates to pixel coordinates using the WCS.
Parameters:
ra (float): Right ascension in degrees.
dec (float): Declination in degrees.
Returns:
Tuple[float, float]: Pixel coordinates (x, y) in the image.
"""
return SkyCoord(ra, dec, unit="deg").to_pixel(self.wcs)
[docs]
def pixels_to_skycoord(self, pixels: np.ndarray):
"""
Convert pixel coordinates to sky coordinates using the WCS.
Parameters:
pixels (np.ndarray): Array of pixel coordinates.
Returns:
Array of sky coordinates corresponding to the input pixels.
"""
return self.wcs.pixel_to_world_values(pixels)
[docs]
def find_gaia_match(self):
"""
Find the closest Gaia star match for each detected star in the image.
Calculates the distance between each detected star and all Gaia stars
in pixel space, then returns the closest match and distance for each.
Returns:
Tuple[np.ndarray, np.ndarray]:
- Matched Gaia star positions in pixel coordinates
- Distances to the closest matches in pixels
"""
squared_distances = np.sum(
(
self.stars_in_image[:, np.newaxis]
- self.gaia_stars_in_image[np.newaxis, :, :]
)
** 2,
axis=-1,
)
match_index = np.argmin(squared_distances, axis=1)
distance = np.sqrt(np.min(squared_distances, axis=1))
return self.gaia_stars_in_image[match_index], distance
[docs]
def number_of_matched_stars(self, pixel_threshold: float = 10):
"""
Count the number of stars matched within a pixel threshold.
Parameters:
pixel_threshold (float, optional): Maximum distance in pixels to consider
a match valid. Defaults to 10.
Returns:
int: Number of detected stars that have a Gaia match within the threshold.
"""
distance_to_closest_star = self.find_gaia_match()[1]
return np.sum(distance_to_closest_star < pixel_threshold)
[docs]
def plot(self, ax=None, transpose=False, **kwargs):
"""
Plot detected stars and their Gaia matches for visualization.
Parameters:
ax (matplotlib.axes.Axes, optional): Axes to plot on. If None, creates new figure.
transpose (bool, optional): If True, transpose x and y coordinates for plotting.
Defaults to False.
**kwargs: Additional keyword arguments passed to scatter plot.
"""
if ax is None:
fig, ax = plt.subplots()
default_dict = {
"facecolors": "none",
"edgecolors": "r",
"linewidth": 2,
"marker": "o",
}
gaia_stars = self.find_gaia_match()
gaia_stars_in_pixel_range = gaia_stars[0][gaia_stars[1] < 20]
dim = (1, 0) if transpose else (0, 1)
ax.scatter(
self.stars_in_image[::, dim[0]],
self.stars_in_image[:, dim[1]],
label="Detected stars",
s=80,
**(default_dict | kwargs),
)
ax.scatter(
gaia_stars_in_pixel_range[:, dim[0]],
gaia_stars_in_pixel_range[:, dim[1]],
label="Matched Gaia stars",
s=160,
**(default_dict | {"edgecolors": "dodgerblue"} | kwargs),
)
ax.scatter(
self.gaia_stars_in_image[:, dim[0]],
self.gaia_stars_in_image[:, dim[1]],
label="All Gaia stars",
s=40,
**(default_dict | {"edgecolors": "yellow"} | kwargs),
)
ax.legend()
def _read_fits_file(filepath: str | Path):
"""
Read image data and header from a FITS file.
Parameters:
filepath (str or Path): Path to the FITS file.
Returns:
Tuple[np.ndarray, Header]: Image data and header.
"""
with fits.open(filepath) as hdul:
header = hdul[0].header
image = hdul[0].data
return image, header
def _extract_plate_scale_and_dateobs(header):
"""
Extract plate scale and observation date from FITS header.
Parameters:
header (fits.Header): FITS header containing metadata.
Returns:
Tuple[float, datetime]: Plate scale in degrees per pixel and observation date.
"""
dateobs = pd.to_datetime(header["DATE-OBS"])
# get units of FOCALLEN through comments if available
focallen_comment = header.comments["FOCALLEN"].lower()
if "mm" in focallen_comment or "millimeter" in focallen_comment:
focallen_unit = 1e-3
else:
focallen_unit = 1.0 # default to m
plate_scale = optics.plate_scale(
header["XPIXSZ"] * 1e-6, header["FOCALLEN"] * focallen_unit
) # deg/pixel
return plate_scale, dateobs
def _map_filter_band_to_gaia_tmass(filter_band: Optional[str]) -> Optional[str]:
"""
Map a given filter band to the corresponding Gaia or 2MASS band.
Parameters:
filter_band (str, optional): The filter band used for observation.
Returns:
str: Corresponding Gaia or 2MASS band, or defaults to G if not found.
"""
gaia_tmass_filter_mappings = {
"G": ["r", "clear", "C"],
"BP": ["u", "g"],
"RP": ["i", "z", "I+z", "Exo"],
"J": ["J", "Y", "YJ", "zYJ"],
"H": ["H"],
"KS": ["K", "Ks"],
}
if filter_band is None:
return "G"
for gaia_band, bands in gaia_tmass_filter_mappings.items():
if filter_band.lower().strip("'") in [b.lower() for b in bands]:
return gaia_band
return "G"
def _get_gaia_star_coordinates(
ra,
dec,
image_clean,
dateobs,
plate_scale,
filter_band=None,
fov_scale=1.1,
limit=24,
use_local_db: bool = False,
logger: Optional[logging.Logger] = None,
):
"""
Get Gaia star coordinates for a given field of view.
Parameters:
ra (float): Right ascension of field center in degrees.
dec (float): Declination of field center in degrees.
image_clean (np.ndarray): The cleaned image data.
dateobs (datetime): Observation date.
plate_scale (float): Plate scale in degrees per pixel.
filter_band (str, optional): Filter band used for observation.
fov_scale (float, optional): Factor to scale field of view. Defaults to 1.1.
limit (int, optional): Maximum number of stars to query. Defaults to 24.
use_local_db (bool, optional): Whether to use the local Gaia database. Defaults to False.
logger (logging.Logger, optional): Logger for progress messages.
Defaults to the module logger.
Returns:
np.ndarray: Array of Gaia star coordinates.
"""
logger = logging.getLogger(__name__) if logger is None else logger
# Get fov from image shape and plate scale
fov = optics.fov_from_plate_scale(image_clean.shape, plate_scale) * fov_scale
gaia_tmass_filter = _map_filter_band_to_gaia_tmass(filter_band)
# Try online query first (if local DB doesn't exist, this is the only option)
has_local_db = Config().gaia_db.is_file()
if use_local_db is False:
try:
table = cabaret.GaiaQuery.query(
center=(ra, dec),
radius=np.max(fov) / 2,
filter_bands=gaia_tmass_filter,
limit=limit,
timeout=60,
)
table_filt = cabaret.GaiaQuery._apply_proper_motion(table, dateobs).copy()
return np.array(
[table_filt["ra"].value.data, table_filt["dec"].value.data]
).T
except Exception as e:
if not has_local_db:
raise Exception(
"Failed to query online Gaia database and no local fallback available"
) from e
logger.warning(
f"Online Gaia query failed with error: {e}. Falling back to local database."
)
# Fall back to local database
use_tmass = gaia_tmass_filter in ["J", "H", "KS"]
logger.debug(
f"Using {'2MASS J' if use_tmass else 'Gaia G'} filter band for local query"
)
return local_gaia_db_query(
center=(ra, dec), fov=fov, limit=limit, dateobs=dateobs, tmass=use_tmass
)
def _verify_offset_within_fov(
pointing_correction: PointingCorrection,
plate_scale: float,
image_shape: Tuple[int, int],
):
"""
Verify that the pointing offset is within the field of view.
Parameters:
pointing_correction (PointingCorrection): The pointing correction object.
plate_scale (float): Plate scale in degrees per pixel.
image_shape (Tuple[int, int]): Shape of the image.
Raises:
Exception: If the offset is larger than the field of view.
"""
# Check that the offset is not larger than the fov
# Get fov from image shape and plate scale
fov = optics.fov_from_plate_scale(image_shape, plate_scale)
if pointing_correction.angular_separation > max(fov):
raise Exception("Pointing error is larger than the field of view")
def _verify_plate_solve(
image_star_mapping: ImageStarMapping,
plate_scale: float,
scale_tolerance: float = 0.05,
max_axis_ratio: float = 1.02,
max_skew_deg: float = 1.0,
match_radius_arcsec: float = 3.5,
min_matched_stars: int = 4,
):
"""
Verify that the WCS from a plate solve is physically plausible.
A wrong solution can match a few stars by chance, but it seldom also has
the correct plate scale and square pixels. The solve is accepted only if
all of these checks pass:
- Each axis of the fitted plate scale is within ``scale_tolerance`` of
``plate_scale``.
- The pixels are square: the ratio of the two axis scales is below
``max_axis_ratio`` and the axes are within ``max_skew_deg`` of 90 degrees.
- At least ``min_matched_stars`` detected stars have a Gaia star within
``match_radius_arcsec``.
Parameters:
image_star_mapping (ImageStarMapping): The image-star mapping object.
plate_scale (float): Expected plate scale in degrees per pixel, for
example from FOCALLEN and XPIXSZ.
scale_tolerance (float, optional): Allowed fractional difference
between the fitted and expected plate scale. Defaults to 0.05.
max_axis_ratio (float, optional): Largest allowed ratio between the two
axis scales. Defaults to 1.02.
max_skew_deg (float, optional): Largest allowed difference in degrees
between the angle of the two axes and 90 degrees. Defaults to 1.0.
match_radius_arcsec (float, optional): Largest distance in arcseconds
between a detected star and its Gaia match. Defaults to 3.5, which
is about 10 pixels at 0.35"/pixel.
min_matched_stars (int, optional): Smallest number of detected stars
that must have a Gaia match. Defaults to 4.
Returns:
str: A one-line summary of the measured values, for the log.
Raises:
PlateSolveError: If one of the checks fails.
"""
wcs = image_star_mapping.wcs
fitted_scales = proj_plane_pixel_scales(wcs) # deg/pixel, one per axis
scale_errors = fitted_scales / plate_scale - 1
if np.any(np.abs(scale_errors) > scale_tolerance):
raise PlateSolveError(
f'fitted plate scale {fitted_scales[0] * 3600:.3f}"/px x '
f'{fitted_scales[1] * 3600:.3f}"/px is not within '
f'{scale_tolerance:.0%} of the expected {plate_scale * 3600:.3f}"/px'
)
axis_ratio = fitted_scales.max() / fitted_scales.min()
cd = wcs.pixel_scale_matrix
cos_angle = np.dot(cd[:, 0], cd[:, 1]) / np.prod(fitted_scales)
skew_deg = np.degrees(np.arcsin(min(abs(cos_angle), 1.0)))
if axis_ratio >= max_axis_ratio or skew_deg >= max_skew_deg:
raise PlateSolveError(
f"pixels are not square: axis ratio {axis_ratio:.3f} "
f"(limit {max_axis_ratio}), skew {skew_deg:.2f} deg "
f"(limit {max_skew_deg} deg)"
)
match_radius_px = match_radius_arcsec / 3600 / plate_scale
number_of_matched_stars = image_star_mapping.number_of_matched_stars(
pixel_threshold=match_radius_px
)
if number_of_matched_stars < min_matched_stars:
raise PlateSolveError(
f"only {number_of_matched_stars} of "
f"{len(image_star_mapping.stars_in_image)} detected stars have a "
f'Gaia star within {match_radius_arcsec}" ({match_radius_px:.1f} px). '
f"Required: {min_matched_stars}."
)
return (
f"plate scale error {scale_errors[0]:+.1%} x {scale_errors[1]:+.1%}, "
f"axis ratio {axis_ratio:.3f}, skew {skew_deg:.2f} deg, "
f"{number_of_matched_stars}/{len(image_star_mapping.stars_in_image)} "
f'stars matched within {match_radius_arcsec}" ({match_radius_px:.1f} px)'
)
[docs]
def local_gaia_db_query(
center: Union[Tuple[float, float], SkyCoord],
fov: Union[float, Quantity],
limit: int = 1000,
tmass: bool = False,
dateobs: Optional[datetime] = None,
) -> np.ndarray:
"""
Query local Gaia database to retrieve RA-DEC coordinates of stars within a given field-of-view.
Retrieves star positions from a local Gaia database within a rectangular region
centered on the specified coordinates. Optionally applies proper motion corrections
for a specific observation date and sorts by magnitude (Gaia G-band or 2MASS J-band).
Parameters:
center (tuple or SkyCoord): The sky coordinates of the center of the FOV.
If a tuple is given, it should contain the RA and DEC in degrees.
fov (float or Quantity): The field-of-view size in degrees. If a float is given,
it is assumed to be in degrees. Can be a single value (square FOV) or
tuple (RA_fov, Dec_fov).
limit (int, optional): The maximum number of sources to retrieve from the
database. Defaults to 1000.
tmass (bool, optional): Whether to sort by 2MASS J magnitudes instead of
Gaia G magnitudes. Defaults to False.
dateobs (datetime, optional): The date of the observation. If given, the
proper motions of the sources will be applied to update positions
from J2015.5 to the observation date. Defaults to None.
Returns:
np.ndarray: An array of shape (n, 2) containing the RA-DEC coordinates
of the retrieved sources in degrees, sorted by magnitude (brightest first).
"""
if isinstance(center, SkyCoord):
ra = center.ra.deg
dec = center.dec.deg
else:
ra, dec = center
logger.debug(f"Querying local Gaia DB around RA={ra}, DEC={dec} with FOV={fov}")
if not isinstance(fov, u.Quantity):
fov = fov * u.deg
logger.debug(f"Converted FOV to Quantity: {fov}")
if fov.ndim == 0:
ra_fov = dec_fov = fov.to(u.deg).value
elif fov.ndim == 1:
ra_fov, dec_fov = fov.to(u.deg).value
else:
ra_fov = fov[0].to(u.deg).value
dec_fov = fov[1].to(u.deg).value
# Dec bounds are straightforward
min_dec = max(dec - dec_fov / 2, -90.0)
max_dec = min(dec + dec_fov / 2, 90.0)
# For RA, account for spherical geometry
# At declination `dec`, the RA angular size scales as 1/cos(dec)
cos_dec = np.cos(np.radians(dec))
# Check if we're near a pole (within 5 degrees)
if abs(dec) > 85:
# Near poles, query all RA values
min_ra = 0.0
max_ra = 360.0
logger.debug("Near pole - querying all RA values")
else:
# Adjust RA FOV for spherical geometry
ra_fov_adjusted = ra_fov / cos_dec if cos_dec > 0.01 else 360.0
min_ra = ra - ra_fov_adjusted / 2
max_ra = ra + ra_fov_adjusted / 2
# Handle RA wraparound at 0/360
if min_ra < 0:
min_ra += 360
if max_ra > 360:
max_ra -= 360
# Handle RA wraparound case
crosses_zero = min_ra > max_ra
table = local_db_query(
Config().gaia_db, min_dec, max_dec, min_ra, max_ra, crosses_zero=crosses_zero
)
if tmass:
table = table.sort_values(by=["j_m"]).reset_index(drop=True)
else:
table = table.sort_values(by=["phot_g_mean_mag"]).reset_index(drop=True)
table = table.map(lambda x: np.nan if x == "" else x)
table.dropna(inplace=True)
# Limit number of stars
table = table[0:limit]
# Add proper motion to ra and dec
if dateobs is not None:
dateobs = dateobs.year + (dateobs.timetuple().tm_yday - 1) / 365.25 # type: ignore
years = dateobs - 2015.5 # type: ignore
table["ra"] += years * table["pmra"] / 1000 / 3600
table["dec"] += years * table["pmdec"] / 1000 / 3600
return np.array([table["ra"].values, table["dec"].values]).T
[docs]
def local_db_query(
db: Union[str, Path],
min_dec: float,
max_dec: float,
min_ra: float,
max_ra: float,
crosses_zero: bool = False,
) -> pd.DataFrame:
"""Query astronomical database for objects within coordinate bounds.
Performs federated query across sharded SQLite database tables to retrieve
astronomical catalog data within specified declination and right ascension ranges.
Args:
db (Union[str, Path]): Path to the SQLite database file.
min_dec (float): Minimum declination in degrees.
max_dec (float): Maximum declination in degrees.
min_ra (float): Minimum right ascension in degrees.
max_ra (float): Maximum right ascension in degrees.
crosses_zero (bool, optional): Whether the RA range crosses the 0-degree line.
Defaults to False.
Returns:
pd.DataFrame: Combined results from all relevant database shards.
"""
conn = sqlite3.connect(db)
# Determine the relevant shard(s) based on declination
arr = np.arange(np.floor(min_dec), np.ceil(max_dec) + 1, 1)
relevant_shard_ids = set()
for i in range(len(arr) - 1):
shard_id = f"{arr[i]:.0f}_{arr[i + 1]:.0f}"
relevant_shard_ids.add(shard_id)
# Execute the federated query across the relevant shard(s)
df_total = pd.DataFrame()
for shard_id in relevant_shard_ids:
shard_table_name = f"{shard_id}"
if crosses_zero:
# Query in two parts: [min_ra, 360] OR [0, max_ra]
q = f"""SELECT * FROM `{shard_table_name}`
WHERE dec BETWEEN {min_dec} AND {max_dec}
AND (ra >= {min_ra} OR ra <= {max_ra})"""
else:
q = f"""SELECT * FROM `{shard_table_name}`
WHERE dec BETWEEN {min_dec} AND {max_dec}
AND ra BETWEEN {min_ra} AND {max_ra}"""
df = pd.read_sql_query(q, conn)
df_total = pd.concat([df, df_total], axis=0)
conn.close()
return df_total