Update to hats healpix math

src/lsdb/dask/merge_catalog_functions.py

@@ -86,7 +86,7 @@ def align_catalogs(left: Catalog, right: Catalog, add_right_margin: bool = True)
         else right.hc_structure.pixel_tree.to_moc()
     )
     if right_added_radius is not None:
-        right_moc_depth_resol = hp.nside2resol(hp.order2nside(right_moc.max_order), arcmin=True) * 60
+        right_moc_depth_resol = hp.order2resol(right_moc.max_order, arcmin=True) * 60
         if right_added_radius < right_moc_depth_resol:
             right_moc = copy_moc(right_moc).add_neighbours()
         else:

src/lsdb/loaders/dataframe/from_dataframe.py

@@ -19,7 +19,7 @@ def from_dataframe(
     drop_empty_siblings: bool = False,
     partition_size: int | None = None,
     threshold: int | None = None,
-    margin_order: int | None = -1,
+    margin_order: int = -1,
     margin_threshold: float | None = 5.0,
     should_generate_moc: bool = True,
     moc_max_order: int = 10,

src/lsdb/loaders/dataframe/margin_catalog_generator.py

@@ -27,7 +27,7 @@ class MarginCatalogGenerator:
     def __init__(
         self,
         catalog: Catalog,
-        margin_order: int | None = -1,
+        margin_order: int = -1,
         margin_threshold: float = 5.0,
         use_pyarrow_types: bool = True,
         **kwargs,
@@ -169,12 +169,10 @@ def _create_margins(self, margin_pairs_df: pd.DataFrame) -> Dict[HealpixPixel, p
             A dictionary mapping each margin pixel to the respective DataFrame.
         """
         margin_pixel_df_map: Dict[HealpixPixel, npd.NestedFrame] = {}
-        self.dataframe["margin_pixel"] = hp.ang2pix(
-            2**self.margin_order,
+        self.dataframe["margin_pixel"] = hp.radec2pix(
+            self.margin_order,
             self.dataframe[self.hc_structure.catalog_info.ra_column].to_numpy(),
             self.dataframe[self.hc_structure.catalog_info.dec_column].to_numpy(),
-            lonlat=True,
-            nest=True,
         )
         constrained_data = self.dataframe.reset_index().merge(margin_pairs_df, on="margin_pixel")
         if len(constrained_data):

tests/lsdb/catalog/test_catalog.py

@@ -10,6 +10,7 @@
 import numpy.testing as npt
 import pandas as pd
 import pytest
+from hats.io.file_io import read_fits_image
 from hats.pixel_math import HealpixPixel, spatial_index_to_healpix
 
 import lsdb
@@ -208,21 +209,13 @@ def test_save_catalog_point_map(small_sky_order1_catalog, tmp_path):
 
     point_map_path = base_catalog_path / "point_map.fits"
     assert hc.io.file_io.does_file_or_directory_exist(point_map_path)
-    map_fits_image = hp.read_map(point_map_path, nest=True, h=True)
-    histogram, header_dict = map_fits_image[0], dict(map_fits_image[1])
+    histogram = read_fits_image(point_map_path)
 
     # The histogram and the sky map histogram match
     assert len(small_sky_order1_catalog) == np.sum(histogram)
     expected_histogram = small_sky_order1_catalog.skymap_histogram(lambda df, _: len(df), order=8)
     npt.assert_array_equal(expected_histogram, histogram)
 
-    # Check the fits file metadata
-    assert header_dict["PIXTYPE"] == "HEALPIX"
-    assert header_dict["ORDERING"] == "NESTED"
-    assert header_dict["INDXSCHM"] == "IMPLICIT"
-    assert header_dict["OBJECT"] == "FULLSKY"
-    assert header_dict["NSIDE"] == 256
-
 
 def test_save_catalog_overwrite(small_sky_catalog, tmp_path):
     base_catalog_path = tmp_path / "small_sky"
@@ -324,7 +317,7 @@ def test_prune_empty_partitions_all_are_removed(small_sky_order1_catalog):
 
 def test_skymap_data(small_sky_order1_catalog):
     def func(df, healpix):
-        return len(df) / hp.nside2pixarea(hp.order2nside(healpix.order), degrees=True)
+        return len(df) / hp.order2pixarea(healpix.order, degrees=True)
 
     skymap = small_sky_order1_catalog.skymap_data(func)
     for pixel in skymap.keys():
@@ -335,7 +328,7 @@ def func(df, healpix):
 
 def test_skymap_data_order(small_sky_order1_catalog):
     def func(df, healpix):
-        return len(df) / hp.nside2pixarea(hp.order2nside(healpix.order), degrees=True)
+        return len(df) / hp.order2pixarea(healpix.order, degrees=True)
 
     order = 3
 
@@ -357,7 +350,7 @@ def func(df, healpix):
 
 def test_skymap_data_wrong_order(small_sky_order1_catalog):
     def func(df, healpix):
-        return len(df) / hp.nside2pixarea(hp.order2nside(healpix.order), degrees=True)
+        return len(df) / hp.order2pixarea(healpix.order, degrees=True)
 
     order = 0
 
@@ -367,7 +360,7 @@ def func(df, healpix):
 
 def test_skymap_histogram(small_sky_order1_catalog):
     def func(df, healpix):
-        return len(df) / hp.nside2pixarea(hp.order2nside(healpix.order), degrees=True)
+        return len(df) / hp.order2pixarea(healpix.order, degrees=True)
 
     pixel_map = small_sky_order1_catalog.skymap_data(func)
     pixel_map = {pixel: value.compute() for pixel, value in pixel_map.items()}
@@ -384,7 +377,7 @@ def func(df, healpix):
 
 def test_skymap_histogram_empty(small_sky_order1_catalog):
     def func(df, healpix):
-        return len(df) / hp.nside2pixarea(hp.order2nside(healpix.order), degrees=True)
+        return len(df) / hp.order2pixarea(healpix.order, degrees=True)
 
     expected_img = np.full(12, 1)
     img = small_sky_order1_catalog.cone_search(0, 0, 1).skymap_histogram(func, default_value=1)
@@ -396,7 +389,7 @@ def test_skymap_histogram_order_default(small_sky_order1_catalog):
     default = -1.0
 
     def func(df, _):
-        return len(df) / hp.nside2pixarea(hp.order2nside(order), degrees=True)
+        return len(df) / hp.order2pixarea(order, degrees=True)
 
     computed_catalog = small_sky_order1_catalog.compute()
     order_3_pixels = spatial_index_to_healpix(computed_catalog.index.to_numpy(), order)
@@ -412,7 +405,7 @@ def test_skymap_histogram_null_values_order_default(small_sky_order1_catalog):
     default = -1.0
 
     def func(df, healpix):
-        density = len(df) / hp.nside2pixarea(hp.order2nside(healpix.order), degrees=True)
+        density = len(df) / hp.order2pixarea(healpix.order, degrees=True)
         return density if healpix.pixel % 2 == 0 else None
 
     pixels = list(small_sky_order1_catalog._ddf_pixel_map.keys())
@@ -438,7 +431,7 @@ def test_skymap_histogram_null_values_order(small_sky_order1_catalog):
     default = -1.0
 
     def func(df, healpix):
-        density = len(df) / hp.nside2pixarea(hp.order2nside(healpix.order), degrees=True)
+        density = len(df) / hp.order2pixarea(healpix.order, degrees=True)
         return density if healpix.pixel % 2 == 0 else None
 
     pixels = list(small_sky_order1_catalog._ddf_pixel_map.keys())
@@ -462,7 +455,7 @@ def test_skymap_histogram_order_empty(small_sky_order1_catalog):
     order = 3
 
     def func(df, healpix):
-        return len(df) / hp.nside2pixarea(hp.order2nside(healpix.order), degrees=True)
+        return len(df) / hp.order2pixarea(healpix.order, degrees=True)
 
     catalog = small_sky_order1_catalog.cone_search(0, 0, 1)
     _, non_empty_partitions = catalog._get_non_empty_partitions()
@@ -477,7 +470,7 @@ def test_skymap_histogram_order_some_partitions_empty(small_sky_order1_catalog):
     order = 3
 
     def func(df, healpix):
-        return len(df) / hp.nside2pixarea(hp.order2nside(healpix.order), degrees=True)
+        return len(df) / hp.order2pixarea(healpix.order, degrees=True)
 
     catalog = small_sky_order1_catalog.query("ra > 350 and dec < -50")
     _, non_empty_partitions = catalog._get_non_empty_partitions()
@@ -502,7 +495,7 @@ def test_skymap_plot(small_sky_order1_catalog, mocker):
     mocker.patch("lsdb.catalog.dataset.healpix_dataset.plot_healpix_map")
 
     def func(df, healpix):
-        return len(df) / hp.nside2pixarea(hp.order2nside(healpix.order), degrees=True)
+        return len(df) / hp.order2pixarea(healpix.order, degrees=True)
 
     small_sky_order1_catalog.skymap(func)
     pixel_map = small_sky_order1_catalog.skymap_data(func)
@@ -600,9 +593,7 @@ def add_col(df, pixel):
     assert isinstance(mapped, Catalog)
     assert "pix" in mapped.columns
     mapcomp = mapped.compute()
-    pix_col = hp.ang2pix(
-        hp.order2nside(1), mapcomp["ra"].to_numpy(), mapcomp["dec"].to_numpy(), lonlat=True, nest=True
-    )
+    pix_col = hp.radec2pix(1, mapcomp["ra"].to_numpy(), mapcomp["dec"].to_numpy())
     assert np.all(mapcomp["pix"] == pix_col)
 
 

tests/lsdb/loaders/dataframe/test_from_dataframe.py

@@ -95,7 +95,7 @@ def test_partitions_on_map_equal_partitions_in_df(small_sky_order1_df, small_sky
         partition_df = catalog._ddf.partitions[partition_index].compute()
         assert isinstance(partition_df, pd.DataFrame)
         for _, row in partition_df.iterrows():
-            ipix = hp.ang2pix(2**hp_pixel.order, row["ra"], row["dec"], nest=True, lonlat=True)
+            ipix = hp.radec2pix(hp_pixel.order, row["ra"], row["dec"])
             assert ipix == hp_pixel.pixel