Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Migrate polygon search to use mocpy utilities #415

Merged
merged 6 commits into from
Nov 14, 2024
Merged
Show file tree
Hide file tree
Changes from 5 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
70 changes: 1 addition & 69 deletions src/hats/catalog/catalog.py
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand All @@ -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.

Expand Down
72 changes: 71 additions & 1 deletion src/hats/catalog/healpix_dataset/healpix_dataset.py
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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
Expand Down Expand Up @@ -100,6 +109,8 @@
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")

Check warning on line 113 in src/hats/catalog/healpix_dataset/healpix_dataset.py

View check run for this annotation

Codecov / codecov/patch

src/hats/catalog/healpix_dataset/healpix_dataset.py#L113

Added line #L113 was not covered by tests
max_order = (
max(self.moc.max_order, self.pixel_tree.get_max_depth())
if self.moc is not None
Expand All @@ -123,6 +134,63 @@
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.

Expand All @@ -132,6 +200,8 @@
Returns:
A new catalog with only the pixels that overlap with the moc. Note that we reset the total_rows
to 0, as updating would require a scan over the new pixel sizes."""
if len(self.pixel_tree) == 0:
raise ValueError("Cannot filter empty catalog")

Check warning on line 204 in src/hats/catalog/healpix_dataset/healpix_dataset.py

View check run for this annotation

Codecov / codecov/patch

src/hats/catalog/healpix_dataset/healpix_dataset.py#L204

Added line #L204 was not covered by tests
filtered_tree = filter_by_moc(self.pixel_tree, moc)
filtered_moc = self.moc.intersection(moc) if self.moc is not None else None
filtered_catalog_info = self.catalog_info.copy_and_update(total_rows=0)
Expand Down
3 changes: 1 addition & 2 deletions src/hats/pixel_math/box_filter.py
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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 = []
Expand Down
7 changes: 0 additions & 7 deletions src/hats/pixel_math/cone_filter.py

This file was deleted.

24 changes: 0 additions & 24 deletions src/hats/pixel_math/polygon_filter.py

This file was deleted.

59 changes: 34 additions & 25 deletions src/hats/pixel_math/validators.py
Original file line number Diff line number Diff line change
Expand Up @@ -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"""
Expand All @@ -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):
Expand Down Expand Up @@ -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
Expand Down
15 changes: 7 additions & 8 deletions tests/hats/catalog/test_catalog.py
Original file line number Diff line number Diff line change
Expand Up @@ -224,14 +224,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):
Expand Down Expand Up @@ -289,6 +285,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):
Expand Down