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 box filtering to use MOCpy #428

Merged
merged 11 commits into from
Nov 22, 2024
2 changes: 1 addition & 1 deletion pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ dependencies = [
"healpy",
"jproperties",
"matplotlib>=3.3,<3.9",
"mocpy",
"mocpy>=0.17.1",
"numba>=0.58",
"numpy<3",
"pandas",
Expand Down
8 changes: 3 additions & 5 deletions src/hats/catalog/healpix_dataset/healpix_dataset.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
98 changes: 10 additions & 88 deletions src/hats/pixel_math/box_filter.py
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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:
Expand All @@ -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))
16 changes: 0 additions & 16 deletions src/hats/pixel_math/healpix_shim.py
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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


Expand Down
24 changes: 11 additions & 13 deletions src/hats/pixel_math/validators.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Loading