From 2193fecbb0694925066aea681c6924726ffc5ddc Mon Sep 17 00:00:00 2001 From: Sandro Campos Date: Mon, 18 Nov 2024 15:06:59 -0500 Subject: [PATCH 1/9] Use mocpy in box filtering --- .../healpix_dataset/healpix_dataset.py | 8 +- src/hats/pixel_math/box_filter.py | 101 ++---------------- src/hats/pixel_math/validators.py | 31 +++--- 3 files changed, 28 insertions(+), 112 deletions(-) diff --git a/src/hats/catalog/healpix_dataset/healpix_dataset.py b/src/hats/catalog/healpix_dataset/healpix_dataset.py index b5971717..177631bf 100644 --- a/src/hats/catalog/healpix_dataset/healpix_dataset.py +++ b/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 + small arc circles and the declination edges follow great arc circles. Args: ra (Tuple[float, float]): Right ascension range, in degrees diff --git a/src/hats/pixel_math/box_filter.py b/src/hats/pixel_math/box_filter.py index ae33ee45..c1419542 100644 --- a/src/hats/pixel_math/box_filter.py +++ b/src/hats/pixel_math/box_filter.py @@ -1,37 +1,29 @@ from __future__ import annotations from collections.abc import Iterable -from typing import List, 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 sides follow + great arc circles, the declination sides follow small arc circles. Args: - ra (Tuple[float, float]): Right ascension range, in [0,360] degrees - dec (Tuple[float, float]): Declination range, in [-90,90] degrees + ra (tuple[float, float]): Right ascension range, in [0,360] degrees + dec (tuple[float, float]): Declination range, in [-90,90] degrees order (int): Maximum order of the moc to generate the box at 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 +36,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( - nside, 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, lonlat=True) - pixels = hp.query_polygon(nside, vertices, inclusive=True, nest=True) - all_polygon_pixels.append(pixels) - return np.unique(np.concatenate(all_polygon_pixels, 0)) diff --git a/src/hats/pixel_math/validators.py b/src/hats/pixel_math/validators.py index 6e60e432..ee0c88e8 100644 --- a/src/hats/pixel_math/validators.py +++ b/src/hats/pixel_math/validators.py @@ -11,7 +11,7 @@ class ValidatorsErrors(str, Enum): """Error messages for the coordinate validators""" - INVALID_DEC = "declination must be in the -90.0 to 90.0 degree range" + INVALID_DEC = "declination must be in the [-90,90[ degree range" INVALID_RADIUS = "cone radius must be positive" INVALID_NUM_VERTICES = "polygon must contain a minimum of 3 vertices" DUPLICATE_VERTICES = "polygon has duplicated vertices" @@ -35,17 +35,17 @@ def validate_radius(radius_arcsec: float): def validate_declination_values(dec: float | List[float]): - """Validates that declination values are in the [-90,90] degree range + """Validates that declination values are in the [-90,90[ degree range Args: dec (float | List[float]): The declination values to be validated Raises: - ValueError: if declination values are not in the [-90,90] degree range + ValueError: if declination values are not in the [-90,90[ degree range """ dec_values = np.array(dec) lower_bound, upper_bound = -90.0, 90.0 - if not np.all((dec_values >= lower_bound) & (dec_values <= upper_bound)): + if not np.all((dec_values >= lower_bound) & (dec_values < upper_bound)): raise ValueError(ValidatorsErrors.INVALID_DEC.value) @@ -98,25 +98,22 @@ 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 unique pair of values, in degrees + - Declination values must be in ascending order and in 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 or ra[0] == ra[1]: + 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) From d36f5bbe6d0647d2c3c5c5294716aa84d92fb3bd Mon Sep 17 00:00:00 2001 From: Sandro Campos Date: Mon, 18 Nov 2024 16:27:40 -0500 Subject: [PATCH 2/9] Update tests --- tests/hats/catalog/test_catalog.py | 155 +++++++---------------------- tests/hats/conftest.py | 17 ++++ 2 files changed, 55 insertions(+), 117 deletions(-) diff --git a/tests/hats/catalog/test_catalog.py b/tests/hats/catalog/test_catalog.py index 988e04a8..29ee3f80 100644 --- a/tests/hats/catalog/test_catalog.py +++ b/tests/hats/catalog/test_catalog.py @@ -1,6 +1,7 @@ """Tests of catalog functionality""" import os +import re import astropy.units as u import numpy as np @@ -15,7 +16,6 @@ from hats.io.file_io import read_fits_image from hats.loaders import read_hats from hats.pixel_math import HealpixPixel -from hats.pixel_math.box_filter import _generate_ra_strip_moc, generate_box_moc from hats.pixel_math.validators import ValidatorsErrors from hats.pixel_tree.pixel_tree import PixelTree @@ -205,9 +205,9 @@ def test_cone_filter_empty(small_sky_order1_catalog): def test_cone_filter_invalid_cone_center(small_sky_order1_catalog): - with pytest.raises(ValueError, match=ValidatorsErrors.INVALID_DEC): + with pytest.raises(ValueError, match=re.escape(ValidatorsErrors.INVALID_DEC)): small_sky_order1_catalog.filter_by_cone(0, -100, 0.1) - with pytest.raises(ValueError, match=ValidatorsErrors.INVALID_DEC): + with pytest.raises(ValueError, match=re.escape(ValidatorsErrors.INVALID_DEC)): small_sky_order1_catalog.filter_by_cone(0, 100, 0.1) with pytest.raises(ValueError, match=ValidatorsErrors.INVALID_RADIUS): small_sky_order1_catalog.filter_by_cone(0, 10, -1) @@ -269,7 +269,7 @@ def test_polygonal_filter_empty(small_sky_order1_catalog): def test_polygonal_filter_invalid_coordinates(small_sky_order1_catalog): # Declination is over 90 degrees polygon_vertices = [(47.1, -100), (64.5, -100), (64.5, 6.27), (47.1, 6.27)] - with pytest.raises(ValueError, match=ValidatorsErrors.INVALID_DEC): + with pytest.raises(ValueError, match=re.escape(ValidatorsErrors.INVALID_DEC)): small_sky_order1_catalog.filter_by_polygon(polygon_vertices) # Right ascension should wrap, it does not throw an error polygon_vertices = [(470.1, 6), (470.5, 6), (64.5, 10.27), (47.1, 10.27)] @@ -297,30 +297,33 @@ def test_polygonal_filter_invalid_polygon(small_sky_order1_catalog): small_sky_order1_catalog.filter_by_polygon(vertices) -def test_box_filter_ra(small_sky_order1_catalog): - ra = (280, 290) - # The catalog pixels are distributed around the [270,0] degree range. - filtered_catalog = small_sky_order1_catalog.filter_by_box(ra=ra) - +def test_box_filter_ra_and_dec(small_sky_order1_catalog): + # The catalog pixels are distributed around the [-90,0] degree range. + filtered_catalog = small_sky_order1_catalog.filter_by_box(ra=(280, 300), dec=(-30, -20)) filtered_pixels = filtered_catalog.get_healpix_pixels() assert len(filtered_pixels) == 2 - assert filtered_pixels == [HealpixPixel(1, 44), HealpixPixel(1, 46)] + assert filtered_pixels == [HealpixPixel(1, 46), HealpixPixel(1, 47)] - assert (1, 44) in filtered_catalog.pixel_tree assert (1, 46) in filtered_catalog.pixel_tree + assert (1, 47) in filtered_catalog.pixel_tree assert len(filtered_catalog.pixel_tree.pixels[1]) == 2 assert filtered_catalog.catalog_info.total_rows == 0 - assert filtered_catalog.moc is not None - box_moc = generate_box_moc(ra=ra, dec=None, order=small_sky_order1_catalog.get_max_coverage_order()) - assert filtered_catalog.moc == box_moc.intersection(small_sky_order1_catalog.moc) + # Check that the previous filter is the same as intersecting the ra and dec filters + filtered_catalog_ra = small_sky_order1_catalog.filter_by_box(ra=(280, 300)) + filtered_catalog_dec = small_sky_order1_catalog.filter_by_box(dec=(-30, -20)) + filtered_catalog_ra_pixels = filtered_catalog_ra.get_healpix_pixels() + filtered_catalog_dec_pixels = filtered_catalog_dec.get_healpix_pixels() + intersected_pixels = [ + pixel for pixel in filtered_catalog_ra_pixels if pixel in filtered_catalog_dec_pixels + ] + assert filtered_pixels == intersected_pixels def test_box_filter_wrapped_ra(small_sky_order1_catalog): # The catalog pixels are distributed around the [270,0] degree range. - filtered_catalog = small_sky_order1_catalog.filter_by_box(ra=(-10, 10)) - + filtered_catalog = small_sky_order1_catalog.filter_by_box(ra=(-10, 10), dec=(-90, 89)) filtered_pixels = filtered_catalog.get_healpix_pixels() assert len(filtered_pixels) == 2 @@ -339,6 +342,7 @@ def test_box_filter_ra_divisions_edge_cases(small_sky_order1_catalog): # that wide RA ranges (larger than 180 degrees) are correctly handled. # The catalog pixels are distributed around the [270,0] degree range. + dec = (-90, 89) def assert_is_subset_of(catalog, catalog_complement): pixels_catalog = catalog.get_healpix_pixels() @@ -346,138 +350,55 @@ def assert_is_subset_of(catalog, catalog_complement): assert len(pixels_catalog) < len(pixels_catalog_complement) assert all(pixel in pixels_catalog_complement for pixel in pixels_catalog) - filtered_catalog = small_sky_order1_catalog.filter_by_box(ra=(0, 180)) - filtered_catalog_complement = small_sky_order1_catalog.filter_by_box(ra=(180, 0)) + filtered_catalog = small_sky_order1_catalog.filter_by_box(ra=(0, 180), dec=dec) + filtered_catalog_complement = small_sky_order1_catalog.filter_by_box(ra=(180, 0), dec=dec) assert_is_subset_of(filtered_catalog, filtered_catalog_complement) - filtered_catalog = small_sky_order1_catalog.filter_by_box(ra=(10, 50)) - filtered_catalog_complement = small_sky_order1_catalog.filter_by_box(ra=(50, 10)) + filtered_catalog = small_sky_order1_catalog.filter_by_box(ra=(10, 50), dec=dec) + filtered_catalog_complement = small_sky_order1_catalog.filter_by_box(ra=(50, 10), dec=dec) assert_is_subset_of(filtered_catalog, filtered_catalog_complement) - filtered_catalog = small_sky_order1_catalog.filter_by_box(ra=(10, 220)) - filtered_catalog_complement = small_sky_order1_catalog.filter_by_box(ra=(220, 10)) + filtered_catalog = small_sky_order1_catalog.filter_by_box(ra=(10, 220), dec=dec) + filtered_catalog_complement = small_sky_order1_catalog.filter_by_box(ra=(220, 10), dec=dec) assert_is_subset_of(filtered_catalog, filtered_catalog_complement) - filtered_catalog = small_sky_order1_catalog.filter_by_box(ra=(350, 200)) - filtered_catalog_complement = small_sky_order1_catalog.filter_by_box(ra=(200, 350)) + filtered_catalog = small_sky_order1_catalog.filter_by_box(ra=(350, 200), dec=dec) + filtered_catalog_complement = small_sky_order1_catalog.filter_by_box(ra=(200, 350), dec=dec) assert_is_subset_of(filtered_catalog, filtered_catalog_complement) - filtered_catalog = small_sky_order1_catalog.filter_by_box(ra=(50, 200)) - filtered_catalog_complement = small_sky_order1_catalog.filter_by_box(ra=(200, 50)) + filtered_catalog = small_sky_order1_catalog.filter_by_box(ra=(50, 200), dec=dec) + filtered_catalog_complement = small_sky_order1_catalog.filter_by_box(ra=(200, 50), dec=dec) assert_is_subset_of(filtered_catalog, filtered_catalog_complement) -def test_box_filter_ra_pixel_tree_generation(): - """This method tests the pixel tree generation for the ra filter""" - # The catalog pixels are distributed around the [270,0] degree range. - moc = _generate_ra_strip_moc(ra_range=(0, 180), order=1) - moc_complement = _generate_ra_strip_moc(ra_range=(180, 0), order=1) - assert len(moc.flatten()) == len(moc_complement.flatten()) - - moc = _generate_ra_strip_moc(ra_range=(10, 50), order=1) - moc_complement = _generate_ra_strip_moc(ra_range=(50, 10), order=1) - assert len(moc.flatten()) < len(moc_complement.flatten()) - - moc = _generate_ra_strip_moc(ra_range=(10, 220), order=1) - moc_complement = _generate_ra_strip_moc(ra_range=(220, 10), order=1) - assert len(moc_complement.flatten()) < len(moc.flatten()) - - moc = _generate_ra_strip_moc(ra_range=(200, 350), order=1) - moc_complement = _generate_ra_strip_moc(ra_range=(350, 200), order=1) - assert len(moc.flatten()) < len(moc_complement.flatten()) - - moc = _generate_ra_strip_moc(ra_range=(200, 50), order=1) - moc_complement = _generate_ra_strip_moc(ra_range=(50, 200), order=1) - assert len(moc_complement.flatten()) < len(moc.flatten()) - - -def test_box_filter_dec(small_sky_order1_catalog): - # The catalog pixels are distributed around the [-90,0] degree range. - filtered_catalog = small_sky_order1_catalog.filter_by_box(dec=(10, 20)) - assert len(filtered_catalog.get_healpix_pixels()) == 0 - assert len(filtered_catalog.pixel_tree) == 0 - assert filtered_catalog.catalog_info.total_rows == 0 - - filtered_catalog_1 = small_sky_order1_catalog.filter_by_box(dec=(-10, 10)) - filtered_pixels_1 = filtered_catalog_1.get_healpix_pixels() - assert filtered_pixels_1 == [HealpixPixel(1, 47)] - assert (1, 47) in filtered_catalog_1.pixel_tree - assert len(filtered_catalog_1.pixel_tree) == 1 - assert filtered_catalog_1.catalog_info.total_rows == 0 - - filtered_catalog_2 = small_sky_order1_catalog.filter_by_box(dec=(-30, -20)) - filtered_pixels_2 = filtered_catalog_2.get_healpix_pixels() - assert filtered_pixels_2 == [HealpixPixel(1, 45), HealpixPixel(1, 46), HealpixPixel(1, 47)] - assert (1, 45) in filtered_catalog_2.pixel_tree - assert (1, 46) in filtered_catalog_2.pixel_tree - assert (1, 47) in filtered_catalog_2.pixel_tree - assert len(filtered_catalog_2.pixel_tree) == 3 - assert filtered_catalog_2.catalog_info.total_rows == 0 - - -def test_box_filter_ra_and_dec(small_sky_order1_catalog): - # The catalog pixels are distributed around the [-90,0] degree range. - filtered_catalog = small_sky_order1_catalog.filter_by_box(ra=(280, 300), dec=(-30, -20)) - filtered_pixels = filtered_catalog.get_healpix_pixels() - - assert len(filtered_pixels) == 2 - assert filtered_pixels == [HealpixPixel(1, 46), HealpixPixel(1, 47)] - - assert (1, 46) in filtered_catalog.pixel_tree - assert (1, 47) in filtered_catalog.pixel_tree - assert len(filtered_catalog.pixel_tree.pixels[1]) == 2 - assert filtered_catalog.catalog_info.total_rows == 0 - - # Check that the previous filter is the same as intersecting the ra and dec filters - filtered_catalog_ra = small_sky_order1_catalog.filter_by_box(ra=(280, 300)) - filtered_catalog_dec = small_sky_order1_catalog.filter_by_box(dec=(-30, -20)) - filtered_catalog_ra_pixels = filtered_catalog_ra.get_healpix_pixels() - filtered_catalog_dec_pixels = filtered_catalog_dec.get_healpix_pixels() - intersected_pixels = [ - pixel for pixel in filtered_catalog_ra_pixels if pixel in filtered_catalog_dec_pixels - ] - assert filtered_pixels == intersected_pixels - - def test_box_filter_empty(small_sky_order1_catalog): - # It is very difficult to get an empty set of HEALPix with ra for this test catalog - # as its pixels are very close to the South Pole (dec of -90 degrees). In order 1, - # they are very large in area, and easily overlap with any ra region. - filtered_catalog = small_sky_order1_catalog.filter_by_box(ra=(0, 10)) - assert len(filtered_catalog.get_healpix_pixels()) == 2 - assert len(filtered_catalog.pixel_tree) == 2 - - filtered_catalog = small_sky_order1_catalog.filter_by_box(dec=(10, 20)) - assert len(filtered_catalog.get_healpix_pixels()) == 0 - assert len(filtered_catalog.pixel_tree) == 0 - filtered_catalog = small_sky_order1_catalog.filter_by_box(ra=(40, 50), dec=(10, 20)) assert len(filtered_catalog.get_healpix_pixels()) == 0 assert len(filtered_catalog.pixel_tree) == 0 def test_box_filter_invalid_args(small_sky_order1_catalog): - # Some declination values are out of the [-90,90] bounds - with pytest.raises(ValueError, match=ValidatorsErrors.INVALID_DEC): + # Some declination values are out of the [-90,90[ bounds + with pytest.raises(ValueError, match=re.escape(ValidatorsErrors.INVALID_DEC)): small_sky_order1_catalog.filter_by_box(ra=(0, 30), dec=(-100, -70)) # Declination values should be in ascending order with pytest.raises(ValueError, match=ValidatorsErrors.INVALID_RADEC_RANGE): - small_sky_order1_catalog.filter_by_box(dec=(0, -10)) + small_sky_order1_catalog.filter_by_box(ra=(0, 30), dec=(0, -10)) # There are ranges are defined with more than 2 values with pytest.raises(ValueError, match=ValidatorsErrors.INVALID_RADEC_RANGE): - small_sky_order1_catalog.filter_by_box(ra=(0, 30), dec=(-30, -40, 10)) + small_sky_order1_catalog.filter_by_box(ra=(0, 30), dec=(-40, -30, 10)) with pytest.raises(ValueError, match=ValidatorsErrors.INVALID_RADEC_RANGE): small_sky_order1_catalog.filter_by_box(ra=(0, 30, 40), dec=(-40, 10)) # The range values coincide (for ra, values are wrapped) with pytest.raises(ValueError, match=ValidatorsErrors.INVALID_RADEC_RANGE): - small_sky_order1_catalog.filter_by_box(ra=(100, 100)) + small_sky_order1_catalog.filter_by_box(ra=(100, 100), dec=(-40, -30)) with pytest.raises(ValueError, match=ValidatorsErrors.INVALID_RADEC_RANGE): - small_sky_order1_catalog.filter_by_box(ra=(0, 360)) + small_sky_order1_catalog.filter_by_box(ra=(0, 360), dec=(-40, -30)) with pytest.raises(ValueError, match=ValidatorsErrors.INVALID_RADEC_RANGE): - small_sky_order1_catalog.filter_by_box(dec=(50, 50)) + small_sky_order1_catalog.filter_by_box(ra=(0, 50), dec=(50, 50)) # No range values were provided with pytest.raises(ValueError, match=ValidatorsErrors.INVALID_RADEC_RANGE): @@ -594,6 +515,6 @@ def test_catalog_len_is_undetermined(small_sky_order1_catalog): vertices = [(300, -50), (300, -55), (272, -55), (272, -50)] len(small_sky_order1_catalog.filter_by_polygon(vertices)) with pytest.raises(ValueError, match="undetermined"): - len(small_sky_order1_catalog.filter_by_box(ra=(280, 300))) + len(small_sky_order1_catalog.filter_by_box(ra=(280, 300), dec=(0, 30))) with pytest.raises(ValueError, match="undetermined"): len(small_sky_order1_catalog.filter_from_pixel_list([HealpixPixel(0, 11)])) diff --git a/tests/hats/conftest.py b/tests/hats/conftest.py index ccbd60a6..a34c06df 100644 --- a/tests/hats/conftest.py +++ b/tests/hats/conftest.py @@ -1,5 +1,8 @@ +import numpy as np import pytest +import hats.pixel_math.healpix_shim as hp +from hats.catalog import Catalog, TableProperties from hats.loaders import read_hats from hats.pixel_math import HealpixPixel @@ -82,3 +85,17 @@ def pixel_list_breadth_first(): HealpixPixel(1, 45), HealpixPixel(1, 46), ] + + +@pytest.fixture +def full_sky_catalog_order5(): + all_pixels_order5 = [HealpixPixel(5, pix) for pix in np.arange(hp.order2npix(5))] + catalog_info = TableProperties( + catalog_name="full_sky_order5", + catalog_type="object", + ra_column="ra", + dec_column="dec", + total_rows=0, + hats_copyright="LINCC Frameworks 2024", + ) + return Catalog(catalog_info, all_pixels_order5) From 5850ebb5b35efb05d1755eacfb2db0634b32e90d Mon Sep 17 00:00:00 2001 From: Sandro Campos Date: Mon, 18 Nov 2024 16:35:56 -0500 Subject: [PATCH 3/9] Remove full sky coverage test catalog --- tests/hats/conftest.py | 17 ----------------- 1 file changed, 17 deletions(-) diff --git a/tests/hats/conftest.py b/tests/hats/conftest.py index a34c06df..ccbd60a6 100644 --- a/tests/hats/conftest.py +++ b/tests/hats/conftest.py @@ -1,8 +1,5 @@ -import numpy as np import pytest -import hats.pixel_math.healpix_shim as hp -from hats.catalog import Catalog, TableProperties from hats.loaders import read_hats from hats.pixel_math import HealpixPixel @@ -85,17 +82,3 @@ def pixel_list_breadth_first(): HealpixPixel(1, 45), HealpixPixel(1, 46), ] - - -@pytest.fixture -def full_sky_catalog_order5(): - all_pixels_order5 = [HealpixPixel(5, pix) for pix in np.arange(hp.order2npix(5))] - catalog_info = TableProperties( - catalog_name="full_sky_order5", - catalog_type="object", - ra_column="ra", - dec_column="dec", - total_rows=0, - hats_copyright="LINCC Frameworks 2024", - ) - return Catalog(catalog_info, all_pixels_order5) From 3218bc1acf81a324bd42e5297d66b45223f08f9a Mon Sep 17 00:00:00 2001 From: Sandro Campos Date: Mon, 18 Nov 2024 16:57:26 -0500 Subject: [PATCH 4/9] Ensure moc is correct in box filter test --- tests/hats/catalog/test_catalog.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/tests/hats/catalog/test_catalog.py b/tests/hats/catalog/test_catalog.py index 29ee3f80..3f06af6d 100644 --- a/tests/hats/catalog/test_catalog.py +++ b/tests/hats/catalog/test_catalog.py @@ -7,6 +7,7 @@ import numpy as np import pyarrow as pa import pytest +from astropy.coordinates import SkyCoord from mocpy import MOC import hats.pixel_math.healpix_shim as hp @@ -297,7 +298,7 @@ def test_polygonal_filter_invalid_polygon(small_sky_order1_catalog): small_sky_order1_catalog.filter_by_polygon(vertices) -def test_box_filter_ra_and_dec(small_sky_order1_catalog): +def test_box_filter(small_sky_order1_catalog): # The catalog pixels are distributed around the [-90,0] degree range. filtered_catalog = small_sky_order1_catalog.filter_by_box(ra=(280, 300), dec=(-30, -20)) filtered_pixels = filtered_catalog.get_healpix_pixels() @@ -311,14 +312,13 @@ def test_box_filter_ra_and_dec(small_sky_order1_catalog): assert filtered_catalog.catalog_info.total_rows == 0 # Check that the previous filter is the same as intersecting the ra and dec filters - filtered_catalog_ra = small_sky_order1_catalog.filter_by_box(ra=(280, 300)) - filtered_catalog_dec = small_sky_order1_catalog.filter_by_box(dec=(-30, -20)) - filtered_catalog_ra_pixels = filtered_catalog_ra.get_healpix_pixels() - filtered_catalog_dec_pixels = filtered_catalog_dec.get_healpix_pixels() - intersected_pixels = [ - pixel for pixel in filtered_catalog_ra_pixels if pixel in filtered_catalog_dec_pixels - ] - assert filtered_pixels == intersected_pixels + assert filtered_catalog.moc is not None + box_moc = MOC.from_zone( + # SkyCoord([bottom_left_corner, upper_right_corner]) + SkyCoord([[280, -30], [300, -20]], unit="deg"), + max_depth=small_sky_order1_catalog.get_max_coverage_order(), + ) + assert filtered_catalog.moc == box_moc.intersection(small_sky_order1_catalog.moc) def test_box_filter_wrapped_ra(small_sky_order1_catalog): From 155c26580b06baf58091af96ee4126cec05aa44e Mon Sep 17 00:00:00 2001 From: Sandro Campos Date: Mon, 18 Nov 2024 17:20:54 -0500 Subject: [PATCH 5/9] Fix confusion in docstring --- src/hats/catalog/healpix_dataset/healpix_dataset.py | 2 +- src/hats/pixel_math/box_filter.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/hats/catalog/healpix_dataset/healpix_dataset.py b/src/hats/catalog/healpix_dataset/healpix_dataset.py index 177631bf..f3d5f027 100644 --- a/src/hats/catalog/healpix_dataset/healpix_dataset.py +++ b/src/hats/catalog/healpix_dataset/healpix_dataset.py @@ -158,7 +158,7 @@ def filter_by_cone(self, ra: float, dec: float, radius_arcsec: float) -> 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 zone, defined by right ascension and declination ranges. The right ascension edges follow - small arc circles and the declination edges follow great arc circles. + great arc circles and the declination edges follow small arc circles. Args: ra (Tuple[float, float]): Right ascension range, in degrees diff --git a/src/hats/pixel_math/box_filter.py b/src/hats/pixel_math/box_filter.py index c1419542..3b6c28b6 100644 --- a/src/hats/pixel_math/box_filter.py +++ b/src/hats/pixel_math/box_filter.py @@ -9,8 +9,8 @@ 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 sides follow - great arc circles, the declination sides follow small arc circles. + 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 From a055795a737c0816c7ddf6f63414c36de0429430 Mon Sep 17 00:00:00 2001 From: Sandro Campos Date: Mon, 18 Nov 2024 17:59:43 -0500 Subject: [PATCH 6/9] Remove unused healpy methods --- src/hats/pixel_math/healpix_shim.py | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/src/hats/pixel_math/healpix_shim.py b/src/hats/pixel_math/healpix_shim.py index bee05e47..34e2b494 100644 --- a/src/hats/pixel_math/healpix_shim.py +++ b/src/hats/pixel_math/healpix_shim.py @@ -36,21 +36,9 @@ def ang2pix(*args, **kwargs): return hp.ang2pix(*args, **kwargs) -def ring2nest(*args, **kwargs): - return hp.ring2nest(*args, **kwargs) - - ## Query -def query_strip(*args, **kwargs): - return hp.query_strip(*args, **kwargs) - - -def query_polygon(*args, **kwargs): - return hp.query_polygon(*args, **kwargs) - - def get_all_neighbours(*args, **kwargs): return hp.get_all_neighbours(*args, **kwargs) From efa8587050ccc418f82347e96fd721ccc9f3b92f Mon Sep 17 00:00:00 2001 From: Sandro Campos Date: Tue, 19 Nov 2024 16:59:13 -0500 Subject: [PATCH 7/9] Improve error message --- src/hats/pixel_math/box_filter.py | 2 +- src/hats/pixel_math/validators.py | 2 +- tests/hats/catalog/test_catalog.py | 10 ++++------ 3 files changed, 6 insertions(+), 8 deletions(-) diff --git a/src/hats/pixel_math/box_filter.py b/src/hats/pixel_math/box_filter.py index 3b6c28b6..c9b33cd2 100644 --- a/src/hats/pixel_math/box_filter.py +++ b/src/hats/pixel_math/box_filter.py @@ -14,7 +14,7 @@ def generate_box_moc(ra: tuple[float, float], dec: tuple[float, float], order: i Args: ra (tuple[float, float]): Right ascension range, in [0,360] degrees - dec (tuple[float, float]): Declination range, in [-90,90] degrees + dec (tuple[float, float]): Declination range, in [-90,90[ degrees order (int): Maximum order of the moc to generate the box at Returns: diff --git a/src/hats/pixel_math/validators.py b/src/hats/pixel_math/validators.py index 56f567fc..0672baf0 100644 --- a/src/hats/pixel_math/validators.py +++ b/src/hats/pixel_math/validators.py @@ -11,7 +11,7 @@ class ValidatorsErrors(str, Enum): """Error messages for the coordinate validators""" - INVALID_DEC = "declination must be in the [-90,90[ degree range" + INVALID_DEC = "declination must be in the -90 <= dec < 90 degree range" INVALID_RADIUS = "cone radius must be positive" INVALID_NUM_VERTICES = "polygon must contain a minimum of 3 vertices" DUPLICATE_VERTICES = "polygon has duplicated vertices" diff --git a/tests/hats/catalog/test_catalog.py b/tests/hats/catalog/test_catalog.py index a7258c60..110d7154 100644 --- a/tests/hats/catalog/test_catalog.py +++ b/tests/hats/catalog/test_catalog.py @@ -1,7 +1,5 @@ """Tests of catalog functionality""" - import os -import re import astropy.units as u import numpy as np @@ -206,9 +204,9 @@ def test_cone_filter_empty(small_sky_order1_catalog): def test_cone_filter_invalid_cone_center(small_sky_order1_catalog): - with pytest.raises(ValueError, match=re.escape(ValidatorsErrors.INVALID_DEC)): + with pytest.raises(ValueError, match=ValidatorsErrors.INVALID_DEC): small_sky_order1_catalog.filter_by_cone(0, -100, 0.1) - with pytest.raises(ValueError, match=re.escape(ValidatorsErrors.INVALID_DEC)): + with pytest.raises(ValueError, match=ValidatorsErrors.INVALID_DEC): small_sky_order1_catalog.filter_by_cone(0, 100, 0.1) with pytest.raises(ValueError, match=ValidatorsErrors.INVALID_RADIUS): small_sky_order1_catalog.filter_by_cone(0, 10, -1) @@ -270,7 +268,7 @@ def test_polygonal_filter_empty(small_sky_order1_catalog): def test_polygonal_filter_invalid_coordinates(small_sky_order1_catalog): # Declination is over 90 degrees polygon_vertices = [(47.1, -100), (64.5, -100), (64.5, 6.27), (47.1, 6.27)] - with pytest.raises(ValueError, match=re.escape(ValidatorsErrors.INVALID_DEC)): + with pytest.raises(ValueError, match=ValidatorsErrors.INVALID_DEC): small_sky_order1_catalog.filter_by_polygon(polygon_vertices) # Right ascension should wrap, it does not throw an error polygon_vertices = [(470.1, 6), (470.5, 6), (64.5, 10.27), (47.1, 10.27)] @@ -379,7 +377,7 @@ def test_box_filter_empty(small_sky_order1_catalog): def test_box_filter_invalid_args(small_sky_order1_catalog): # Some declination values are out of the [-90,90[ bounds - with pytest.raises(ValueError, match=re.escape(ValidatorsErrors.INVALID_DEC)): + with pytest.raises(ValueError, match=ValidatorsErrors.INVALID_DEC): small_sky_order1_catalog.filter_by_box(ra=(0, 30), dec=(-100, -70)) # Declination values should be in ascending order From dfe0c2d188f722d79ab0420ece0b1331dacb096f Mon Sep 17 00:00:00 2001 From: Sandro Campos Date: Fri, 22 Nov 2024 14:40:41 -0500 Subject: [PATCH 8/9] Pin mocpy --- pyproject.toml | 2 +- src/hats/pixel_math/validators.py | 2 +- tests/hats/catalog/test_catalog.py | 20 +++++++++++++++----- 3 files changed, 17 insertions(+), 7 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 79f3c456..7502b8b7 100644 --- a/pyproject.toml +++ b/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", diff --git a/src/hats/pixel_math/validators.py b/src/hats/pixel_math/validators.py index 0672baf0..88d43c51 100644 --- a/src/hats/pixel_math/validators.py +++ b/src/hats/pixel_math/validators.py @@ -110,7 +110,7 @@ def validate_box(ra: Tuple[float, float], dec: Tuple[float, float]): dec (Tuple[float, float]): Declination range, in degrees """ invalid_range = False - if ra is None or len(ra) != 2 or ra[0] == ra[1]: + 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 diff --git a/tests/hats/catalog/test_catalog.py b/tests/hats/catalog/test_catalog.py index 110d7154..9c749bf9 100644 --- a/tests/hats/catalog/test_catalog.py +++ b/tests/hats/catalog/test_catalog.py @@ -1,4 +1,5 @@ """Tests of catalog functionality""" + import os import astropy.units as u @@ -333,6 +334,19 @@ def test_box_filter_wrapped_ra(small_sky_order1_catalog): assert filtered_catalog.catalog_info.total_rows == 0 +def test_box_filter_ra_boundary(small_sky_order1_catalog): + dec = (-30, 0) + filtered_catalog = small_sky_order1_catalog.filter_by_box(ra=(0, 0), dec=dec) + filtered_pixels = filtered_catalog.get_healpix_pixels() + + assert len(filtered_pixels) == 3 + assert filtered_pixels == [HealpixPixel(1, 45), HealpixPixel(1, 46), HealpixPixel(1, 47)] + + for ra_range in [(0, 360), (360, 0)]: + catalog = small_sky_order1_catalog.filter_by_box(ra=ra_range, dec=dec) + assert catalog.get_healpix_pixels() == filtered_catalog.get_healpix_pixels() + + def test_box_filter_ra_divisions_edge_cases(small_sky_order1_catalog): # In this test we generate RA bands and their complements and compare the amount of # pixels from the catalog after filtering. We construct these complement regions in @@ -390,11 +404,7 @@ def test_box_filter_invalid_args(small_sky_order1_catalog): with pytest.raises(ValueError, match=ValidatorsErrors.INVALID_RADEC_RANGE): small_sky_order1_catalog.filter_by_box(ra=(0, 30, 40), dec=(-40, 10)) - # The range values coincide (for ra, values are wrapped) - with pytest.raises(ValueError, match=ValidatorsErrors.INVALID_RADEC_RANGE): - small_sky_order1_catalog.filter_by_box(ra=(100, 100), dec=(-40, -30)) - with pytest.raises(ValueError, match=ValidatorsErrors.INVALID_RADEC_RANGE): - small_sky_order1_catalog.filter_by_box(ra=(0, 360), dec=(-40, -30)) + # The declination values coincide with pytest.raises(ValueError, match=ValidatorsErrors.INVALID_RADEC_RANGE): small_sky_order1_catalog.filter_by_box(ra=(0, 50), dec=(50, 50)) From 6c36c628023d7ef40ea3e7f236bed296a57c169d Mon Sep 17 00:00:00 2001 From: Sandro Campos Date: Fri, 22 Nov 2024 15:13:08 -0500 Subject: [PATCH 9/9] Test edge cases with dec of 90 deg --- src/hats/pixel_math/box_filter.py | 7 ++++--- src/hats/pixel_math/validators.py | 15 ++++++++------- tests/hats/catalog/test_catalog.py | 6 ++---- 3 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/hats/pixel_math/box_filter.py b/src/hats/pixel_math/box_filter.py index c9b33cd2..ad2d6a13 100644 --- a/src/hats/pixel_math/box_filter.py +++ b/src/hats/pixel_math/box_filter.py @@ -1,20 +1,21 @@ from __future__ import annotations from collections.abc import Iterable +from typing import Tuple import numpy as np from astropy.coordinates import SkyCoord from mocpy import MOC -def generate_box_moc(ra: tuple[float, float], dec: tuple[float, float], order: int) -> MOC: +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 - dec (tuple[float, float]): Declination range, in [-90,90[ degrees + ra (Tuple[float, float]): Right ascension range, in [0,360] degrees + dec (Tuple[float, float]): Declination range, in [-90,90] degrees order (int): Maximum order of the moc to generate the box at Returns: diff --git a/src/hats/pixel_math/validators.py b/src/hats/pixel_math/validators.py index 88d43c51..e81507ca 100644 --- a/src/hats/pixel_math/validators.py +++ b/src/hats/pixel_math/validators.py @@ -11,7 +11,7 @@ class ValidatorsErrors(str, Enum): """Error messages for the coordinate validators""" - INVALID_DEC = "declination must be in the -90 <= dec < 90 degree range" + INVALID_DEC = "declination must be in the -90.0 to 90.0 degree range" INVALID_RADIUS = "cone radius must be positive" INVALID_NUM_VERTICES = "polygon must contain a minimum of 3 vertices" DUPLICATE_VERTICES = "polygon has duplicated vertices" @@ -35,17 +35,17 @@ def validate_radius(radius_arcsec: float): def validate_declination_values(dec: float | List[float]): - """Validates that declination values are in the [-90,90[ degree range + """Validates that declination values are in the [-90,90] degree range Args: dec (float | List[float]): The declination values to be validated Raises: - ValueError: if declination values are not in the [-90,90[ degree range + ValueError: if declination values are not in the [-90,90] degree range """ dec_values = np.array(dec) lower_bound, upper_bound = -90.0, 90.0 - if not np.all((dec_values >= lower_bound) & (dec_values < upper_bound)): + if not np.all((dec_values >= lower_bound) & (dec_values <= upper_bound)): raise ValueError(ValidatorsErrors.INVALID_DEC.value) @@ -101,9 +101,10 @@ def check_polygon_is_valid(vertices: np.ndarray): def validate_box(ra: Tuple[float, float], dec: Tuple[float, float]): """Checks if ra and dec values are valid for the box search. - - Both ranges for ra or dec must have been provided - - Ranges must be defined by a unique pair of values, in degrees - - Declination values must be in ascending order and 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 diff --git a/tests/hats/catalog/test_catalog.py b/tests/hats/catalog/test_catalog.py index 9c749bf9..be96cdcc 100644 --- a/tests/hats/catalog/test_catalog.py +++ b/tests/hats/catalog/test_catalog.py @@ -322,7 +322,7 @@ def test_box_filter(small_sky_order1_catalog): def test_box_filter_wrapped_ra(small_sky_order1_catalog): # The catalog pixels are distributed around the [270,0] degree range. - filtered_catalog = small_sky_order1_catalog.filter_by_box(ra=(-10, 10), dec=(-90, 89)) + filtered_catalog = small_sky_order1_catalog.filter_by_box(ra=(-10, 10), dec=(-90, 90)) filtered_pixels = filtered_catalog.get_healpix_pixels() assert len(filtered_pixels) == 2 @@ -352,9 +352,7 @@ def test_box_filter_ra_divisions_edge_cases(small_sky_order1_catalog): # pixels from the catalog after filtering. We construct these complement regions in # a way that allows us to capture more pixels of the catalog. This is useful to test # that wide RA ranges (larger than 180 degrees) are correctly handled. - - # The catalog pixels are distributed around the [270,0] degree range. - dec = (-90, 89) + dec = (-90, 90) def assert_is_subset_of(catalog, catalog_complement): pixels_catalog = catalog.get_healpix_pixels()