Migrate box filtering to use MOCpy (#428)

* Use mocpy in box filtering

* Update tests

* Ensure moc is correct in box filter test

* Remove unused healpy methods

* Pin mocpy

pyproject.toml

@@ -25,7 +25,7 @@ dependencies = [
     "healpy",
     "jproperties",
     "matplotlib>=3.3,<3.9",
-    "mocpy",
+    "mocpy>=0.17.1",
     "numba>=0.58",
     "numpy<3", 
     "pandas",

src/hats/catalog/healpix_dataset/healpix_dataset.py

@@ -155,12 +155,10 @@ def filter_by_cone(self, ra: float, dec: float, radius_arcsec: float) -> Self:
         )
         return self.filter_by_moc(cone_moc)
 
-    def filter_by_box(
-        self, ra: Tuple[float, float] | None = None, dec: Tuple[float, float] | None = None
-    ) -> Self:
+    def filter_by_box(self, ra: Tuple[float, float], dec: Tuple[float, float]) -> Self:
         """Filter the pixels in the catalog to only include the pixels that overlap with a
-        right ascension or declination range. In case both ranges are provided, filtering
-        is performed using a polygon.
+        zone, defined by right ascension and declination ranges. The right ascension edges follow
+        great arc circles and the declination edges follow small arc circles.
 
         Args:
             ra (Tuple[float, float]): Right ascension range, in degrees

src/hats/pixel_math/box_filter.py

@@ -1,16 +1,17 @@
 from __future__ import annotations
 
 from collections.abc import Iterable
-from typing import List, Tuple
+from typing import Tuple
 
 import numpy as np
+from astropy.coordinates import SkyCoord
 from mocpy import MOC
 
-import hats.pixel_math.healpix_shim as hp
 
-
-def generate_box_moc(ra: Tuple[float, float] | None, dec: Tuple[float, float] | None, order: int) -> MOC:
-    """Generates a MOC object that covers the specified box area
+def generate_box_moc(ra: Tuple[float, float], dec: Tuple[float, float], order: int) -> MOC:
+    """Generates a MOC object that covers the specified box. A box is delimited
+    by right ascension and declination ranges. The right ascension edges follow
+    great arc circles and the declination edges follow small arc circles.
 
     Args:
         ra (Tuple[float, float]): Right ascension range, in [0,360] degrees
@@ -20,18 +21,10 @@ def generate_box_moc(ra: Tuple[float, float] | None, dec: Tuple[float, float] |
     Returns:
         a MOC object that covers the specified box
     """
-    filter_moc = None
-    ra_search_moc, dec_search_moc = None, None
-
-    if ra is not None:
-        ra_search_moc = _generate_ra_strip_moc(ra, order)
-        filter_moc = ra_search_moc
-    if dec is not None:
-        dec_search_moc = _generate_dec_strip_moc(dec, order)
-        filter_moc = dec_search_moc
-    if ra_search_moc is not None and dec_search_moc is not None:
-        filter_moc = ra_search_moc.intersection(dec_search_moc)
-    return filter_moc
+    bottom_left_corner = [ra[0], min(dec)]
+    upper_right_corner = [ra[1], max(dec)]
+    box_coords = SkyCoord([bottom_left_corner, upper_right_corner], unit="deg")
+    return MOC.from_zone(box_coords, max_depth=order)
 
 
 def wrap_ra_angles(ra: np.ndarray | Iterable | int | float) -> np.ndarray:
@@ -44,74 +37,3 @@ def wrap_ra_angles(ra: np.ndarray | Iterable | int | float) -> np.ndarray:
         A numpy array of right ascension values, wrapped to the [0,360] degree range.
     """
     return np.asarray(ra, dtype=float) % 360
-
-
-def _generate_ra_strip_moc(ra_range: Tuple[float, float], order: int) -> MOC:
-    """Generates a pixel_tree filled with leaf nodes at a given order that overlap with the ra region."""
-    # Subdivide polygons (if needed) in two smaller polygons of at most 180 degrees
-    division_ra = _get_division_ra(ra_range)
-
-    if division_ra is not None:
-        pixels_in_range = _get_pixels_for_subpolygons(
-            [
-                # North Pole
-                [(0, 90), (ra_range[0], 0), (division_ra, 0)],
-                [(0, 90), (division_ra, 0), (ra_range[1], 0)],
-                # South Pole
-                [(ra_range[0], 0), (0, -90), (division_ra, 0)],
-                [(division_ra, 0), (0, -90), (ra_range[1], 0)],
-            ],
-            order,
-        )
-    else:
-        # We assume the polygon of smallest area by default
-        # e.g. for ra_ranges of [10,50], [200,10] or [350,10]
-        pixels_in_range = _get_pixels_for_subpolygons(
-            [
-                # North Pole
-                [(0, 90), (ra_range[0], 0), (ra_range[1], 0)],
-                # South Pole
-                [(ra_range[0], 0), (0, -90), (ra_range[1], 0)],
-            ],
-            order,
-        )
-
-    orders = np.full(pixels_in_range.shape, fill_value=order)
-    return MOC.from_healpix_cells(ipix=pixels_in_range, depth=orders, max_depth=order)
-
-
-def _generate_dec_strip_moc(dec_range: Tuple[float, float], order: int) -> MOC:
-    """Generates a pixel_tree filled with leaf nodes at a given order that overlap with the dec region."""
-    nside = hp.order2nside(order)
-    # Convert declination values to colatitudes, in radians, and revert their order
-    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(
-        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)
-
-
-def _get_division_ra(ra_range: Tuple[float, float]) -> float | None:
-    """Gets the division angle for the right ascension region. This is crucial for the edge
-    cases where we need to split up polygons wider than 180 degrees into smaller polygons."""
-    division_ra = None
-    if ra_range[0] > ra_range[1] and 360 - ra_range[0] + ra_range[1] >= 180:
-        # e.g. edge case of [350, 170] or [180, 0], we want the wider polygon
-        division_ra = (ra_range[1] - 360 + ra_range[0]) / 2
-    elif ra_range[0] < ra_range[1] and ra_range[1] - ra_range[0] >= 180:
-        # e.g. edge case of [10, 200] or [0, 180], we want the wider polygon
-        division_ra = ra_range[0] + (ra_range[1] - ra_range[0]) / 2
-    return division_ra
-
-
-def _get_pixels_for_subpolygons(polygons: List[List[tuple[float, float]]], order: int) -> np.ndarray:
-    """Gets the unique pixels for a set of sub-polygons."""
-    nside = hp.order2nside(order)
-    all_polygon_pixels = []
-    for vertices in polygons:
-        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

@@ -4,7 +4,6 @@
 
 import astropy.units as u
 import cdshealpix
-import healpy as hp
 import numpy as np
 from astropy.coordinates import Latitude, Longitude, SkyCoord
 
@@ -68,21 +67,6 @@ def radec2pix(order: int, ra: float, dec: float) -> np.ndarray[np.int64]:
     return cdshealpix.lonlat_to_healpix(ra, dec, order).astype(np.int64)
 
 
-def ring2nest(order: int, ipix: int) -> int:
-    return cdshealpix.from_ring(ipix, order)
-
-
-## Query
-
-
-def query_strip(*args, **kwargs):
-    return hp.query_strip(*args, **kwargs)
-
-
-def query_polygon(*args, **kwargs):
-    return hp.query_polygon(*args, **kwargs)
-
-
 ## Coordinate conversion
 
 

src/hats/pixel_math/validators.py

@@ -98,25 +98,23 @@ def check_polygon_is_valid(vertices: np.ndarray):
             raise ValueError(ValidatorsErrors.INVALID_CONCAVE_SHAPE.value)
 
 
-def validate_box(ra: Tuple[float, float] | None, dec: Tuple[float, float] | None):
+def validate_box(ra: Tuple[float, float], dec: Tuple[float, float]):
     """Checks if ra and dec values are valid for the box search.
 
-    - At least one range of ra or dec must have been provided
-    - Ranges must be pairs of non-duplicate minimum and maximum values, in degrees
-    - Declination values, if existing, must be in ascending order
-    - Declination values, if existing, must be in the [-90,90] degree range
+    - Both ranges for ra or dec must have been provided.
+    - Ranges must be defined by a pair of values, in degrees.
+    - Declination values must be unique, provided in ascending order, and
+    belong to the [-90,90] degree range.
 
     Args:
         ra (Tuple[float, float]): Right ascension range, in degrees
         dec (Tuple[float, float]): Declination range, in degrees
     """
     invalid_range = False
-    if ra is not None:
-        if len(ra) != 2 or len(ra) != len(set(ra)):
-            invalid_range = True
-    if dec is not None:
-        if len(dec) != 2 or dec[0] >= dec[1]:
-            invalid_range = True
-        validate_declination_values(list(dec))
-    if (ra is None and dec is None) or invalid_range:
+    if ra is None or len(ra) != 2:
+        invalid_range = True
+    elif dec is None or len(dec) != 2 or dec[0] >= dec[1]:
+        invalid_range = True
+    if invalid_range:
         raise ValueError(ValidatorsErrors.INVALID_RADEC_RANGE.value)
+    validate_declination_values(dec)