Skip to content

Commit

Permalink
Merge 98bc347 into 17ec4ac
Browse files Browse the repository at this point in the history
  • Loading branch information
camposandro authored Nov 8, 2024
2 parents 17ec4ac + 98bc347 commit d8baa5f
Show file tree
Hide file tree
Showing 9 changed files with 131 additions and 136 deletions.
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
70 changes: 69 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 __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
Expand All @@ -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.
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
2 changes: 2 additions & 0 deletions src/hats/pixel_tree/moc_filter.py
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down
22 changes: 14 additions & 8 deletions tests/hats/catalog/test_catalog.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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):
Expand Down
Loading

0 comments on commit d8baa5f

Please sign in to comment.