import astropy.units as u
import numpy as np
from astropy.coordinates import SkyCoord
from cdshealpix import skycoord_to_healpix
from mocpy import MOC
from scipy.spatial import KDTree
[docs]
def ra_dec_to_cartesian(ra, dec):
"""
Batch convert right ascension and declination to Cartesian coordinates.
We use this custom function over Astropy's built-in conversion for performance reasons.
Because we know the data is in degrees and we only need Cartesian coordinates, we can skip
object creation and units. The results are roughly 50x faster.
Parameters
----------
ra: float or numpy.ndarray
Right ascension in degrees.
dec: float or numpy.ndarray
Declination in degrees.
Returns
-------
x: float or numpy.ndarray
X coordinate.
y: float or numpy.ndarray
Y coordinate.
z: float or numpy.ndarray
Z coordinate.
"""
# Check the bounds of the inputs and handle wrapping in RA.
ra = np.asarray(ra) % 360.0
dec = np.asarray(dec)
if np.any(dec < -90.0) or np.any(dec > 90.0):
raise ValueError("Declination values must be in the range [-90, 90] degrees.")
ra_rad = np.radians(ra)
dec_rad = np.radians(dec)
x = np.cos(dec_rad) * np.cos(ra_rad)
y = np.cos(dec_rad) * np.sin(ra_rad)
z = np.sin(dec_rad)
return x, y, z
[docs]
def dedup_coords(ra, dec, threshold=1e-5):
"""
Remove duplicate coordinates within a specified threshold.
Parameters
----------
ra: numpy.ndarray
Array of right ascension values in degrees.
dec: numpy.ndarray
Array of declination values in degrees.
threshold: float
Minimum separation in degrees to consider two points as distinct.
Returns
-------
unique_ra: numpy.ndarray
Array of unique right ascension values.
unique_dec: numpy.ndarray
Array of unique declination values.
unique_indices: numpy.ndarray
Indices of the unique coordinates in the original arrays.
"""
ra = np.asarray(ra)
dec = np.asarray(dec)
if np.size(ra) != np.size(dec):
raise ValueError("RA and Dec arrays must have the same length.")
# Create a KD-tree for efficient nearest neighbor search.
x, y, z = ra_dec_to_cartesian(ra, dec)
cart_coords = np.array([x, y, z]).T
kd_tree = KDTree(cart_coords)
# Do a range search with the same points to find all neighbors.
adjusted_radius = 2.0 * np.sin(0.5 * np.radians(threshold))
close_points = kd_tree.query_ball_point(cart_coords, adjusted_radius)
# Find unique coordinates. We keep the first occurrence of each unique point.
# Note there will always be at least one match (the point itself).
unique_indices = []
for idx, matches in enumerate(close_points):
if len(matches) == 1 or idx == np.min(matches):
unique_indices.append(idx)
unique_indices = np.array(unique_indices)
unique_ra = ra[unique_indices]
unique_dec = dec[unique_indices]
return unique_ra, unique_dec, unique_indices
[docs]
def build_moc_from_coords(ra, dec, depth=8):
"""
Build a MOC at a given depth from RA and Dec coordinates, such that each
point is covered by the MOC.
Parameters
----------
ra: numpy.ndarray
Array of right ascension values in degrees.
dec: numpy.ndarray
Array of declination values in degrees.
depth: int
Maximum depth of the MOC. Higher values give finer resolution but larger MOC size.
Returns
-------
moc: MOC
MOC object representing the sky coverage of the input coordinates.
"""
if len(ra) != len(dec):
raise ValueError("RA and Dec arrays must have the same length.")
coords = SkyCoord(ra=ra * u.deg, dec=dec * u.deg, frame="icrs")
hpix = skycoord_to_healpix(coords, depth)
unique_hpix = np.unique(hpix)
return MOC.from_healpix_cells(unique_hpix, depth=depth, max_depth=depth)