Migrate polygon search to use mocpy utilities (#415)

* Refactor search

* Improve docstrings

* Add convexity check

* Add more unit tests

src/hats/catalog/catalog.py

@@ -2,22 +2,10 @@
 
 from __future__ import annotations
 
-from typing import List, Tuple
+from typing import List
 
-import numpy as np
-
-import hats.pixel_math.healpix_shim as hp
 from hats.catalog.healpix_dataset.healpix_dataset import HealpixDataset
 from hats.pixel_math import HealpixPixel
-from hats.pixel_math.box_filter import generate_box_moc, wrap_ra_angles
-from hats.pixel_math.cone_filter import generate_cone_moc
-from hats.pixel_math.polygon_filter import CartesianCoordinates, SphericalCoordinates, generate_polygon_moc
-from hats.pixel_math.validators import (
-    validate_box_search,
-    validate_declination_values,
-    validate_polygon,
-    validate_radius,
-)
 from hats.pixel_tree.negative_tree import compute_negative_tree_pixels
 
 
@@ -29,62 +17,6 @@ class Catalog(HealpixDataset):
     `Norder=/Dir=/Npix=.parquet`
     """
 
-    def filter_by_cone(self, ra: float, dec: float, radius_arcsec: float) -> Catalog:
-        """Filter the pixels in the catalog to only include the pixels that overlap with a cone
-
-        Args:
-            ra (float): Right Ascension of the center of the cone in degrees
-            dec (float): Declination of the center of the cone in degrees
-            radius_arcsec (float): Radius of the cone in arcseconds
-
-        Returns:
-            A new catalog with only the pixels that overlap with the specified cone
-        """
-        validate_radius(radius_arcsec)
-        validate_declination_values(dec)
-        return self.filter_by_moc(generate_cone_moc(ra, dec, radius_arcsec, self.get_max_coverage_order()))
-
-    def filter_by_box(
-        self, ra: Tuple[float, float] | None = None, dec: Tuple[float, float] | None = None
-    ) -> Catalog:
-        """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.
-
-        Args:
-            ra (Tuple[float, float]): Right ascension range, in degrees
-            dec (Tuple[float, float]): Declination range, in degrees
-
-        Returns:
-            A new catalog with only the pixels that overlap with the specified region
-        """
-        ra = tuple(wrap_ra_angles(ra)) if ra else None
-        validate_box_search(ra, dec)
-        return self.filter_by_moc(generate_box_moc(ra, dec, self.get_max_coverage_order()))
-
-    def filter_by_polygon(self, vertices: List[SphericalCoordinates] | List[CartesianCoordinates]) -> Catalog:
-        """Filter the pixels in the catalog to only include the pixels that overlap
-        with a polygonal sky region.
-
-        Args:
-            vertices (List[SphericalCoordinates] | List[CartesianCoordinates]): The vertices
-                of the polygon to filter points with, in lists of (ra,dec) or (x,y,z) points
-                on the unit sphere.
-
-        Returns:
-            A new catalog with only the pixels that overlap with the specified polygon.
-        """
-        if all(len(vertex) == 2 for vertex in vertices):
-            ra, dec = np.array(vertices).T
-            validate_declination_values(dec)
-            # Get the coordinates vector on the unit sphere if we were provided
-            # with polygon spherical coordinates of ra and dec
-            cart_vertices = hp.ang2vec(ra, dec, lonlat=True)
-        else:
-            cart_vertices = vertices
-        validate_polygon(cart_vertices)
-        return self.filter_by_moc(generate_polygon_moc(cart_vertices, self.get_max_coverage_order()))
-
     def generate_negative_tree_pixels(self) -> List[HealpixPixel]:
         """Get the leaf nodes at each healpix order that have zero catalog data.
 

src/hats/catalog/healpix_dataset/healpix_dataset.py

@@ -1,11 +1,13 @@
 from __future__ import annotations
 
 from pathlib import Path
-from typing import List, Union
+from typing import List, Tuple, Union
 
+import astropy.units as u
 import numpy as np
 import pandas as pd
 import pyarrow as pa
+from astropy.coordinates import SkyCoord
 from mocpy import MOC
 from typing_extensions import Self
 from upath import UPath
@@ -16,6 +18,13 @@
 from hats.inspection import plot_pixels
 from hats.inspection.visualize_catalog import plot_moc
 from hats.pixel_math import HealpixPixel
+from hats.pixel_math.box_filter import generate_box_moc, wrap_ra_angles
+from hats.pixel_math.validators import (
+    validate_box,
+    validate_declination_values,
+    validate_polygon,
+    validate_radius,
+)
 from hats.pixel_tree import PixelAlignment, PixelAlignmentType
 from hats.pixel_tree.moc_filter import filter_by_moc
 from hats.pixel_tree.pixel_alignment import align_with_mocs
@@ -100,6 +109,8 @@ def __len__(self):
     def get_max_coverage_order(self) -> int:
         """Gets the maximum HEALPix order for which the coverage of the catalog is known from the pixel
         tree and moc if it exists"""
+        if len(self.pixel_tree) == 0:
+            raise ValueError("Cannot get max_order of empty catalog")
         max_order = (
             max(self.moc.max_order, self.pixel_tree.get_max_depth())
             if self.moc is not None
@@ -123,6 +134,63 @@ def filter_from_pixel_list(self, pixels: List[HealpixPixel]) -> Self:
         moc = MOC.from_healpix_cells(ipix=pixel_inds, depth=orders, max_depth=max_order)
         return self.filter_by_moc(moc)
 
+    def filter_by_cone(self, ra: float, dec: float, radius_arcsec: float) -> Self:
+        """Filter the pixels in the catalog to only include the pixels that overlap with a cone
+
+        Args:
+            ra (float): Right ascension of the center of the cone, in degrees
+            dec (float): Declination of the center of the cone, in degrees
+            radius_arcsec (float): Radius of the cone, in arcseconds
+
+        Returns:
+            A new catalog with only the pixels that overlap with the specified cone
+        """
+        validate_radius(radius_arcsec)
+        validate_declination_values(dec)
+        cone_moc = MOC.from_cone(
+            lon=ra * u.deg,
+            lat=dec * u.deg,
+            radius=radius_arcsec * u.arcsec,
+            max_depth=self.get_max_coverage_order(),
+        )
+        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:
+        """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.
+
+        Args:
+            ra (Tuple[float, float]): Right ascension range, in degrees
+            dec (Tuple[float, float]): Declination range, in degrees
+
+        Returns:
+            A new catalog with only the pixels that overlap with the specified region
+        """
+        ra = tuple(wrap_ra_angles(ra)) if ra else None
+        validate_box(ra, dec)
+        box_moc = generate_box_moc(ra, dec, self.get_max_coverage_order())
+        return self.filter_by_moc(box_moc)
+
+    def filter_by_polygon(self, vertices: list[tuple[float, float]]) -> Self:
+        """Filter the pixels in the catalog to only include the pixels that overlap
+        with a polygonal sky region.
+
+        Args:
+            vertices (list[tuple[float,float]]): The list of vertice coordinates for
+                the polygon, (ra, dec), in degrees.
+
+        Returns:
+            A new catalog with only the pixels that overlap with the specified polygon.
+        """
+        validate_polygon(vertices)
+        polygon_moc = MOC.from_polygon_skycoord(
+            SkyCoord(vertices, unit="deg"), max_depth=self.get_max_coverage_order()
+        )
+        return self.filter_by_moc(polygon_moc)
+
     def filter_by_moc(self, moc: MOC) -> Self:
         """Filter the pixels in the catalog to only include the pixels that overlap with the moc provided.
 

src/hats/pixel_math/box_filter.py

@@ -7,7 +7,6 @@
 from mocpy import MOC
 
 import hats.pixel_math.healpix_shim as hp
-from hats.pixel_math.polygon_filter import SphericalCoordinates
 
 
 def generate_box_moc(ra: Tuple[float, float] | None, dec: Tuple[float, float] | None, order: int) -> MOC:
@@ -107,7 +106,7 @@ def _get_division_ra(ra_range: Tuple[float, float]) -> float | None:
     return division_ra
 
 
-def _get_pixels_for_subpolygons(polygons: List[List[SphericalCoordinates]], order: int) -> np.ndarray:
+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 = []

src/hats/pixel_math/validators.py

@@ -5,6 +5,8 @@
 
 import numpy as np
 
+import hats.pixel_math.healpix_shim as hp
+
 
 class ValidatorsErrors(str, Enum):
     """Error messages for the coordinate validators"""
@@ -15,6 +17,8 @@ class ValidatorsErrors(str, Enum):
     DUPLICATE_VERTICES = "polygon has duplicated vertices"
     DEGENERATE_POLYGON = "polygon is degenerate"
     INVALID_RADEC_RANGE = "invalid ra or dec range"
+    INVALID_COORDS_SHAPE = "invalid coordinates shape"
+    INVALID_CONCAVE_SHAPE = "polygon must be convex"
 
 
 def validate_radius(radius_arcsec: float):
@@ -45,51 +49,56 @@ def validate_declination_values(dec: float | List[float]):
         raise ValueError(ValidatorsErrors.INVALID_DEC.value)
 
 
-def validate_polygon(vertices: np.ndarray):
+def validate_polygon(vertices: list[tuple[float, float]]):
     """Checks if the polygon contain a minimum of three vertices, that they are
     unique and that the polygon does not fall on a great circle.
 
     Args:
-        vertices (np.ndarray): The polygon vertices, in cartesian coordinates
+        vertices (list[tuple[float,float]]): The list of vertice coordinates for
+            the polygon, (ra, dec), in degrees.
 
     Raises:
         ValueError: exception if the polygon is invalid.
     """
+    vertices = np.array(vertices)
+    if vertices.shape[1] != 2:
+        raise ValueError(ValidatorsErrors.INVALID_COORDS_SHAPE.value)
+    _, dec = vertices.T
+    validate_declination_values(dec)
     if len(vertices) < 3:
         raise ValueError(ValidatorsErrors.INVALID_NUM_VERTICES.value)
     if len(vertices) != len(np.unique(vertices, axis=0)):
         raise ValueError(ValidatorsErrors.DUPLICATE_VERTICES.value)
-    if is_polygon_degenerate(vertices):
-        raise ValueError(ValidatorsErrors.DEGENERATE_POLYGON.value)
+    check_polygon_is_valid(vertices)
 
 
-def is_polygon_degenerate(vertices: np.ndarray) -> bool:
-    """Checks if all the vertices of the polygon are contained in a same plane.
-    If the plane intersects the center of the sphere, the polygon is degenerate.
+def check_polygon_is_valid(vertices: np.ndarray):
+    """Check if the polygon has no degenerate corners and it is convex.
+
+    Based on HEALpy's `queryPolygonInternal` implementation:
+    https://github.com/cds-astro/cds.moc/blob/master/src/healpix/essentials/HealpixBase.java.
 
     Args:
         vertices (np.ndarray): The polygon vertices, in cartesian coordinates
 
     Returns:
-        A boolean, which is True if the polygon is degenerate, i.e. if it falls
-        on a great circle, False otherwise.
+        True if polygon is valid, False otherwise.
     """
-    # Calculate the normal vector of the plane using three of the vertices
-    normal_vector = np.cross(vertices[1] - vertices[0], vertices[2] - vertices[0])
-
-    # Check if the other vertices lie on the same plane
-    for vertex in vertices[3:]:
-        dot_product = np.dot(normal_vector, vertex - vertices[0])
-        if not np.isclose(dot_product, 0):
-            return False
-
-    # Check if the plane intersects the sphere's center. If it does,
-    # the polygon is degenerate and therefore, invalid.
-    center_distance = np.dot(normal_vector, vertices[0])
-    return bool(np.isclose(center_distance, 0))
-
-
-def validate_box_search(ra: Tuple[float, float] | None, dec: Tuple[float, float] | None):
+    vertices_xyz = hp.ang2vec(*vertices.T, lonlat=True)
+    n_vertices = len(vertices_xyz)
+    flip = 0
+    for i in range(n_vertices):
+        normal = np.cross(vertices_xyz[i], vertices_xyz[(i + 1) % n_vertices])
+        hnd = normal.dot(vertices_xyz[(i + 2) % n_vertices])
+        if np.isclose(hnd, 0, atol=1e-10):
+            raise ValueError(ValidatorsErrors.DEGENERATE_POLYGON.value)
+        if i == 0:
+            flip = -1 if hnd < 0 else 1
+        elif flip * hnd <= 0:
+            raise ValueError(ValidatorsErrors.INVALID_CONCAVE_SHAPE.value)
+
+
+def validate_box(ra: Tuple[float, float] | None, dec: Tuple[float, float] | None):
     """Checks if ra and dec values are valid for the box search.
 
     - At least one range of ra or dec must have been provided

src/hats/pixel_tree/moc_filter.py

@@ -18,6 +18,8 @@ def filter_by_moc(
     Returns:
         A new PixelTree object with only the pixels from the input tree that overlap with the moc.
     """
+    if len(tree) == 0:
+        return tree
     moc_ranges = moc.to_depth29_ranges
     # Convert tree intervals to order 29 to match moc intervals
     tree_29_ranges = tree.tree << (2 * (29 - tree.tree_order))

tests/hats/catalog/test_catalog.py

@@ -10,6 +10,7 @@
 
 import hats.pixel_math.healpix_shim as hp
 from hats.catalog import Catalog, PartitionInfo, TableProperties
+from hats.catalog.healpix_dataset.healpix_dataset import HealpixDataset
 from hats.io import paths
 from hats.io.file_io import read_fits_image
 from hats.loaders import read_hats
@@ -147,6 +148,12 @@ def test_max_coverage_order(small_sky_order1_catalog):
     )
 
 
+def test_max_coverage_order_empty_catalog(catalog_info):
+    empty_catalog = HealpixDataset(catalog_info, PixelTree.from_healpix([]))
+    with pytest.raises(ValueError, match="empty catalog"):
+        empty_catalog.get_max_coverage_order()
+
+
 def test_cone_filter(small_sky_order1_catalog):
     ra = 315
     dec = -66.443
@@ -224,14 +231,10 @@ def test_polygonal_filter(small_sky_order1_catalog):
     assert filtered_catalog.moc == polygon_moc.intersection(small_sky_order1_catalog.moc)
 
 
-def test_polygonal_filter_with_cartesian_coordinates(small_sky_order1_catalog):
-    sky_vertices = [(282, -58), (282, -55), (272, -55), (272, -58)]
-    cartesian_vertices = hp.ang2vec(*np.array(sky_vertices).T, lonlat=True)
-    filtered_catalog_1 = small_sky_order1_catalog.filter_by_polygon(sky_vertices)
-    filtered_catalog_2 = small_sky_order1_catalog.filter_by_polygon(cartesian_vertices)
-    assert filtered_catalog_1.get_healpix_pixels() == filtered_catalog_2.get_healpix_pixels()
-    assert (1, 46) in filtered_catalog_1.pixel_tree
-    assert (1, 46) in filtered_catalog_2.pixel_tree
+def test_polygonal_filter_invalid_coordinate_shape(small_sky_order1_catalog):
+    with pytest.raises(ValueError, match="coordinates shape"):
+        vertices = [(282, -58, 1), (282, -55, 2), (272, -55, 3)]
+        small_sky_order1_catalog.filter_by_polygon(vertices)
 
 
 def test_polygonal_filter_big(small_sky_order1_catalog):
@@ -289,6 +292,9 @@ def test_polygonal_filter_invalid_polygon(small_sky_order1_catalog):
     with pytest.raises(ValueError, match=ValidatorsErrors.DEGENERATE_POLYGON):
         vertices = [(50.1, 0), (100.1, 0), (150.1, 0), (200.1, 0)]
         small_sky_order1_catalog.filter_by_polygon(vertices)
+    with pytest.raises(ValueError, match=ValidatorsErrors.INVALID_CONCAVE_SHAPE):
+        vertices = [(45, 30), (60, 60), (90, 45), (60, 50)]
+        small_sky_order1_catalog.filter_by_polygon(vertices)
 
 
 def test_box_filter_ra(small_sky_order1_catalog):