Migrate ang2vec

src/hats/pixel_math/box_filter.py

@@ -111,7 +111,7 @@ def _get_pixels_for_subpolygons(polygons: List[List[tuple[float, float]]], order
     nside = hp.order2nside(order)
     all_polygon_pixels = []
     for vertices in polygons:
-        vertices = hp.ang2vec(*np.array(vertices).T, lonlat=True)
+        vertices = hp.ang2vec(*np.array(vertices).T)
         pixels = hp.query_polygon(nside, vertices, inclusive=True, nest=True)
         all_polygon_pixels.append(pixels)
     return np.unique(np.concatenate(all_polygon_pixels, 0))

src/hats/pixel_math/healpix_shim.py

@@ -1,7 +1,9 @@
 from __future__ import annotations
 
+import astropy.units as u
 import healpy as hp
 import numpy as np
+from astropy.coordinates import SkyCoord
 
 # pylint: disable=missing-function-docstring
 
@@ -50,8 +52,10 @@ def query_polygon(*args, **kwargs):
 ## Coordinate conversion
 
 
-def ang2vec(*args, **kwargs):
-    return hp.ang2vec(*args, **kwargs)
+def ang2vec(ra, dec, **kwargs) -> np.ndarray:
+    """Converts ra and dec to cartesian coordinates on the unit sphere"""
+    coords = SkyCoord(ra=ra * u.deg, dec=dec * u.deg, **kwargs).cartesian
+    return np.array([coords.x.value, coords.y.value, coords.z.value]).T
 
 
 ## FITS

src/hats/pixel_math/validators.py

@@ -84,7 +84,7 @@ def check_polygon_is_valid(vertices: np.ndarray):
     Returns:
         True if polygon is valid, False otherwise.
     """
-    vertices_xyz = hp.ang2vec(*vertices.T, lonlat=True)
+    vertices_xyz = hp.ang2vec(*vertices.T)
     n_vertices = len(vertices_xyz)
     flip = 0
     for i in range(n_vertices):

tests/hats/pixel_math/test_healpix_shim.py

@@ -37,3 +37,16 @@ def test_order2mindist():
     assert_allclose(hps.order2mindist(orders), min_distances, rtol=1e-2)
 
     assert_allclose(hps.order2mindist(17), 0.01677, rtol=1e-2)
+
+
+def test_ang2vec():
+    """Tests conversion of ra/dec to unit cartesian vectors"""
+    ra = [230.14467816, 110.40507118, 9.41764689, 245.5553117]
+    dec = [38.78080888, 17.09584081, -28.6352765, 5.41803306]
+    ra_rad = np.radians(ra)
+    dec_rad = np.radians(dec)
+    x = np.cos(ra_rad) * np.cos(dec_rad)
+    y = np.sin(ra_rad) * np.cos(dec_rad)
+    z = np.sin(dec_rad)
+    actual = np.asarray([x, y, z]).T
+    assert_array_equal(actual, hps.ang2vec(ra, dec))