Merge pull request #430 from astronomy-commons/sean/healpy-math

Replace healpix with cdshealpix for pixel math operations

src/hats/catalog/partition_info.py

@@ -274,7 +274,7 @@ def calculate_fractional_coverage(self):
         """Calculate what fraction of the sky is covered by partition tiles."""
         pixel_orders = [p.order for p in self.pixel_list]
         cov_order, cov_count = np.unique(pixel_orders, return_counts=True)
-        area_by_order = [hp.nside2pixarea(hp.order2nside(order), degrees=True) for order in cov_order]
+        area_by_order = [hp.order2pixarea(order, degrees=True) for order in cov_order]
         # 41253 is the number of square degrees in a sphere
         # https://en.wikipedia.org/wiki/Square_degree
         return (area_by_order * cov_count).sum() / (360**2 / np.pi)

src/hats/pixel_math/box_filter.py

@@ -87,7 +87,7 @@ def _generate_dec_strip_moc(dec_range: Tuple[float, float], order: int) -> MOC:
     colat_rad = np.sort(np.radians([90 - dec if dec > 0 else 90 + abs(dec) for dec in dec_range]))
     # Figure out which pixels in nested order are contained in the declination band
     pixels_in_range = hp.ring2nest(
-        nside, hp.query_strip(nside, theta1=colat_rad[0], theta2=colat_rad[1], inclusive=True)
+        order, hp.query_strip(nside, theta1=colat_rad[0], theta2=colat_rad[1], inclusive=True)
     )
     orders = np.full(pixels_in_range.shape, fill_value=order)
     return MOC.from_healpix_cells(ipix=pixels_in_range, depth=orders, max_depth=order)

src/hats/pixel_math/healpix_shim.py

@@ -1,41 +1,75 @@
 from __future__ import annotations
 
+import math
+
 import astropy.units as u
+import cdshealpix
 import healpy as hp
 import numpy as np
-from astropy.coordinates import SkyCoord
+from astropy.coordinates import Latitude, Longitude, SkyCoord
 
 # pylint: disable=missing-function-docstring
 
 ## Arithmetic conversions
 
 
-def npix2order(param):
-    return hp.nside2order(hp.npix2nside(param))
+MAX_HEALPIX_ORDER = 29
+
+
+def is_order_valid(order: int) -> bool:
+    return np.all(0 <= order) and np.all(order <= MAX_HEALPIX_ORDER)
+
+
+def npix2order(npix: int) -> int:
+    if npix <= 0:
+        raise ValueError("Invalid value for npix")
+    order = int(math.log2(npix / 12)) >> 1
+    if not is_order_valid(order) or not 12 * (1 << (2 * order)) == npix:
+        raise ValueError("Invalid value for npix")
+    return order
+
+
+def order2nside(order: int) -> int:
+    if not is_order_valid(order):
+        raise ValueError("Invalid value for order")
+    return 1 << order
+
+
+def order2npix(order: int) -> int:
+    if not is_order_valid(order):
+        raise ValueError("Invalid value for order")
+    return 12 * (1 << (2 * order))
 
 
-def order2nside(param):
-    return hp.order2nside(param)
+def order2resol(order: int, arcmin: bool = False) -> float:
+    resol_rad = np.sqrt(order2pixarea(order))
 
+    if arcmin:
+        return np.rad2deg(resol_rad) * 60
 
-def order2npix(param):
-    return hp.order2npix(param)
+    return resol_rad
 
 
-def nside2resol(*args, **kwargs):
-    return hp.nside2resol(*args, **kwargs)
+def order2pixarea(order: int, degrees: bool = False) -> float:
+    npix = order2npix(order)
+    pix_area_rad = 4 * np.pi / npix
+    if degrees:
+        return pix_area_rad * (180 / np.pi) * (180 / np.pi)
+    return pix_area_rad
 
 
-def nside2pixarea(*args, **kwargs):
-    return hp.nside2pixarea(*args, **kwargs)
+def radec2pix(order: int, ra: float, dec: float) -> int:
+    if not is_order_valid(order):
+        raise ValueError("Invalid value for order")
 
+    ra = Longitude(ra, unit="deg")
+    dec = Latitude(dec, unit="deg")
 
-def ang2pix(*args, **kwargs):
-    return hp.ang2pix(*args, **kwargs)
+    return cdshealpix.lonlat_to_healpix(ra, dec, order)
 
 
-def ring2nest(*args, **kwargs):
-    return hp.ring2nest(*args, **kwargs)
+def ring2nest(order: int, ipix: int) -> int:
+    return cdshealpix.from_ring(ipix, order)
 
 
 ## Query
@@ -176,6 +210,5 @@ def order2mindist(order: np.ndarray | int) -> np.ndarray | float:
     np.ndarray of float or a single scalar float
         The minimum distance between pixels in arcminutes
     """
-    pixel_nside = order2nside(order)
-    pixel_avgsize = nside2resol(pixel_nside, arcmin=True)
+    pixel_avgsize = order2resol(order, arcmin=True)
     return avgsize2mindist(pixel_avgsize)

src/hats/pixel_math/partition_stats.py

@@ -43,12 +43,10 @@ def generate_histogram(
     required_columns = [ra_column, dec_column]
     if not all(x in data.columns for x in required_columns):
         raise ValueError(f"Invalid column names in input: {ra_column}, {dec_column}")
-    mapped_pixels = hp.ang2pix(
-        2**highest_order,
+    mapped_pixels = hp.radec2pix(
+        highest_order,
         data[ra_column].values,
         data[dec_column].values,
-        lonlat=True,
-        nest=True,
     )
     mapped_pixel, count_at_pixel = np.unique(mapped_pixels, return_counts=True)
     histogram_result[mapped_pixel] += count_at_pixel.astype(np.int64)

src/hats/pixel_math/spatial_index.py

@@ -24,7 +24,7 @@ def compute_spatial_index(ra_values: List[float], dec_values: List[float]) -> np
     if len(ra_values) != len(dec_values):
         raise ValueError("ra and dec arrays should have the same length")
 
-    mapped_pixels = hp.ang2pix(2**SPATIAL_INDEX_ORDER, ra_values, dec_values, nest=True, lonlat=True)
+    mapped_pixels = hp.radec2pix(SPATIAL_INDEX_ORDER, ra_values, dec_values)
 
     return mapped_pixels
 

tests/hats/pixel_math/test_healpix_shim.py

@@ -1,4 +1,7 @@
+import cdshealpix
 import numpy as np
+import pytest
+from astropy.coordinates import Latitude, Longitude
 from numpy.testing import assert_allclose, assert_array_equal
 
 from hats.pixel_math import healpix_shim as hps
@@ -19,8 +22,7 @@ def test_mindist2avgsize2mindist():
 def test_order2avgsize2order():
     """Test that avgsize2order is inverse of hps.nside2resol(hps.order2nside, arcmin=True)"""
     order = np.arange(20)
-    nside = hps.order2nside(order)
-    assert_array_equal(hps.avgsize2order(hps.nside2resol(nside, arcmin=True)), order)
+    assert_array_equal(hps.avgsize2order(hps.order2resol(order, arcmin=True)), order)
 
 
 def test_margin2order():
@@ -50,3 +52,109 @@ def test_ang2vec():
     z = np.sin(dec_rad)
     actual = np.asarray([x, y, z]).T
     assert_array_equal(actual, hps.ang2vec(ra, dec))
+
+
+def test_npix2order():
+    orders = [0, 1, 5, 10, 20, 29]
+    npix = [12 * (4**order) for order in orders]
+    test_orders = [hps.npix2order(x) for x in npix]
+    assert test_orders == orders
+
+
+def test_npix2order_invalid():
+    npixs = [-10, 0, 11, 13, 47, 49, 100000, 100000000000000000]
+    for npix in npixs:
+        with pytest.raises(ValueError, match="Invalid"):
+            hps.npix2order(npix)
+
+
+def test_order2nside():
+    orders = [0, 1, 5, 10, 20, 29]
+    expected_nsides = [2**x for x in orders]
+    test_nsides = [hps.order2nside(o) for o in orders]
+    assert test_nsides == expected_nsides
+
+
+def test_order2nside_invalid():
+    orders = [-1000, -1, 30, 4000]
+    for order in orders:
+        with pytest.raises(ValueError, match="Invalid"):
+            hps.order2nside(order)
+
+
+def test_order2npix():
+    orders = [0, 1, 5, 10, 20, 29]
+    npix = [12 * (4**order) for order in orders]
+    test_npix = [hps.order2npix(o) for o in orders]
+    assert test_npix == npix
+
+
+def test_order2npix_invalid():
+    orders = [-1000, -1, 30, 4000]
+    for order in orders:
+        with pytest.raises(ValueError, match="Invalid"):
+            hps.order2npix(order)
+
+
+def test_order2pixarea():
+    orders = [0, 1, 5, 10, 20, 29]
+    npix = [12 * (4**order) for order in orders]
+    pix_area_expected = [4 * np.pi / x for x in npix]
+    pix_area_test = [hps.order2pixarea(order) for order in orders]
+    assert pix_area_test == pix_area_expected
+
+
+def test_order2pixarea_degrees():
+    orders = [0, 1, 5, 10, 20, 29]
+    npix = [12 * (4**order) for order in orders]
+    pix_area_expected = [np.rad2deg(np.rad2deg(4 * np.pi / x)) for x in npix]
+    pix_area_test = [hps.order2pixarea(order, degrees=True) for order in orders]
+    assert pix_area_test == pix_area_expected
+
+
+def test_order2resol():
+    orders = [0, 1, 5, 10, 20, 29]
+    resol_expected = [np.sqrt(hps.order2pixarea(order)) for order in orders]
+    resol_test = [hps.order2resol(order) for order in orders]
+    assert resol_test == resol_expected
+
+
+def test_order2resol_arcmin():
+    orders = [0, 1, 5, 10, 20, 29]
+    resol_expected = [np.rad2deg(np.sqrt(hps.order2pixarea(order))) * 60 for order in orders]
+    resol_test = [hps.order2resol(order, arcmin=True) for order in orders]
+    assert resol_test == resol_expected
+
+
+def test_radec2pix_lonlat():
+    orders = [0, 1, 5, 10, 20, 29]
+    ras = np.arange(-180.0, 180.0, 10.0)
+    decs = np.arange(-90.0, 90.0, 180 // len(ras))
+    for order in orders:
+        expected_pixels = cdshealpix.lonlat_to_healpix(
+            Longitude(ras, unit="deg"), Latitude(decs, unit="deg"), order
+        )
+        pixels = hps.radec2pix(order, ras, decs)
+        assert np.all(pixels == expected_pixels)
+
+
+def test_radec2pix_invalid():
+    orders = [0, 1, 5, 10, 20, 29]
+    invalid_orders = [-1000, -1, 30, 40]
+    ras = np.arange(-4000.0, 1000.0, 100.0)
+    decs = np.arange(-1000.0, 1000.0, 2000.0 // len(ras))
+    for order in invalid_orders:
+        with pytest.raises(ValueError, match="Invalid"):
+            hps.radec2pix(order, ras, decs)
+    for order in orders:
+        with pytest.raises(ValueError, match="angle"):
+            hps.radec2pix(order, ras, decs)
+
+
+def test_ring2nest():
+    orders = [0, 1, 5, 10, 20, 29]
+    for order in orders:
+        ipix = np.arange(0, 12 * (4**order), 12 * (4**order) // 10)
+        expected_ring_ipix = cdshealpix.from_ring(ipix, order)
+        test_ring_ipix = hps.ring2nest(order, ipix)
+        assert np.all(test_ring_ipix == expected_ring_ipix)

tests/hats/pixel_math/test_spatial_index.py

@@ -16,12 +16,10 @@
 def test_single():
     """Single point. Adheres to specification."""
     result = compute_spatial_index([5], [5])
-    expected = hp.ang2pix(
-        2**SPATIAL_INDEX_ORDER,
+    expected = hp.radec2pix(
+        SPATIAL_INDEX_ORDER,
         [5],
         [5],
-        nest=True,
-        lonlat=True,
     )
 
     npt.assert_array_equal(result, expected)
@@ -38,12 +36,10 @@ def test_short_list():
     ra = [5, 1, 5]
     dec = [5, 1, 5]
     result = compute_spatial_index(ra, dec)
-    expected = hp.ang2pix(
-        2**SPATIAL_INDEX_ORDER,
+    expected = hp.radec2pix(
+        SPATIAL_INDEX_ORDER,
         ra,
         dec,
-        nest=True,
-        lonlat=True,
     )
     npt.assert_array_equal(result, expected)
 
@@ -53,12 +49,10 @@ def test_list():
     ra = [5, 5, 5, 1, 5, 5, 5, 1, 5]
     dec = [5, 5, 5, 1, 5, 5, 5, 1, 5]
     result = compute_spatial_index(ra, dec)
-    expected = hp.ang2pix(
-        2**SPATIAL_INDEX_ORDER,
+    expected = hp.radec2pix(
+        SPATIAL_INDEX_ORDER,
         ra,
         dec,
-        nest=True,
-        lonlat=True,
     )
     npt.assert_array_equal(result, expected)