From 03b42dc70de207e66b53327f4e8e125f86191e6c Mon Sep 17 00:00:00 2001 From: Melissa DeLucchi Date: Mon, 18 Nov 2024 10:50:54 -0500 Subject: [PATCH 01/19] Remove margin fine filtering. --- .github/workflows/pre-commit-ci.yml | 2 +- .github/workflows/testing-and-coverage.yml | 2 +- src/hats/pixel_math/__init__.py | 1 - src/hats/pixel_math/healpix_shim.py | 20 --- src/hats/pixel_math/margin_bounding.py | 139 ------------------ tests/hats/pixel_math/test_margin_bounding.py | 109 -------------- tests/hats/pixel_math/test_spatial_index.py | 19 ++- 7 files changed, 11 insertions(+), 281 deletions(-) delete mode 100644 src/hats/pixel_math/margin_bounding.py delete mode 100644 tests/hats/pixel_math/test_margin_bounding.py diff --git a/.github/workflows/pre-commit-ci.yml b/.github/workflows/pre-commit-ci.yml index f31b3ee6..642e7b9d 100644 --- a/.github/workflows/pre-commit-ci.yml +++ b/.github/workflows/pre-commit-ci.yml @@ -7,7 +7,7 @@ on: push: branches: [ main ] pull_request: - branches: [ main ] + branches: [ main, margin ] jobs: pre-commit-ci: diff --git a/.github/workflows/testing-and-coverage.yml b/.github/workflows/testing-and-coverage.yml index 3d3b867f..8438b944 100644 --- a/.github/workflows/testing-and-coverage.yml +++ b/.github/workflows/testing-and-coverage.yml @@ -7,7 +7,7 @@ on: push: branches: [ main ] pull_request: - branches: [ main ] + branches: [ main, margin ] jobs: build: diff --git a/src/hats/pixel_math/__init__.py b/src/hats/pixel_math/__init__.py index 132eaa84..327d83f9 100644 --- a/src/hats/pixel_math/__init__.py +++ b/src/hats/pixel_math/__init__.py @@ -2,7 +2,6 @@ from .healpix_pixel import HealpixPixel from .healpix_pixel_convertor import HealpixInputTypes, get_healpix_pixel -from .margin_bounding import check_margin_bounds from .partition_stats import empty_histogram, generate_alignment, generate_histogram from .pixel_margins import get_margin from .spatial_index import compute_spatial_index, spatial_index_to_healpix diff --git a/src/hats/pixel_math/healpix_shim.py b/src/hats/pixel_math/healpix_shim.py index 766bb5d7..44aed38b 100644 --- a/src/hats/pixel_math/healpix_shim.py +++ b/src/hats/pixel_math/healpix_shim.py @@ -20,18 +20,10 @@ def npix2nside(param): return hp.npix2nside(param) -def nside2npix(param): - return hp.nside2npix(param) - - def order2npix(param): return hp.order2npix(param) -def npix2order(param): - return hp.npix2order(param) - - def nside2resol(*args, **kwargs): return hp.nside2resol(*args, **kwargs) @@ -63,10 +55,6 @@ def query_polygon(*args, **kwargs): return hp.query_polygon(*args, **kwargs) -def boundaries(*args, **kwargs): - return hp.boundaries(*args, **kwargs) - - def get_all_neighbours(*args, **kwargs): return hp.get_all_neighbours(*args, **kwargs) @@ -78,14 +66,6 @@ def ang2vec(*args, **kwargs): return hp.ang2vec(*args, **kwargs) -def pix2ang(*args, **kwargs): - return hp.pix2ang(*args, **kwargs) - - -def vec2dir(*args, **kwargs): - return hp.vec2dir(*args, **kwargs) - - def pix2xyf(*args, **kwargs): return hp.pix2xyf(*args, **kwargs) diff --git a/src/hats/pixel_math/margin_bounding.py b/src/hats/pixel_math/margin_bounding.py deleted file mode 100644 index 8234f740..00000000 --- a/src/hats/pixel_math/margin_bounding.py +++ /dev/null @@ -1,139 +0,0 @@ -"""Utilities to build bounding boxes around healpixels that include a neighor margin.""" - -import numpy as np -from astropy.coordinates import SkyCoord -from numba import njit - -import hats.pixel_math.healpix_shim as hp - - -def check_margin_bounds(r_asc, dec, pixel_order, pixel, margin_threshold, step=100, chunk_size=1000): - # pylint: disable=too-many-locals - """Check a set of points in ra and dec space to see if they fall within the - `margin_threshold` of a given healpixel. - - Args: - r_asc (:obj:`np.array`): one dimmensional array representing the ra of a - given set of points. - dec (:obj:`np.array`): one dimmensional array representing the dec of a - given set of points. - pixel_order (int): the order of the healpixel to be checked against. - pixel (int): the healpixel to be checked against. - margin_threshold (double): the size of the margin cache in arcseconds. - step (int): the amount of samples we'll get from the healpixel boundary - to test the datapoints against. the higher the step the more - granular the point checking and therefore the more accurate, however - using a smaller step value might be helpful for super large datasets - to save compute time if accuracy is less important. - chunk_size (int): the size of the chunk of data points we process on a single - thread at any given time. We only process a certain amount of the data - at once as trying to check all data points at once can lead to - memory issues. - Returns: - :obj:`numpy.array` of boolean values corresponding to the ra and dec - coordinates checked against whether a given point is within any of the - provided polygons. - """ - num_points = len(r_asc) - if num_points != len(dec): - raise ValueError( - f"length of r_asc ({num_points}) is different \ - from length of dec ({len(dec)})" - ) - - if num_points == 0: - return np.array([]) - - bounds = hp.vec2dir(hp.boundaries(2**pixel_order, pixel, step=step, nest=True), lonlat=True) - margin_check = np.full((num_points,), False) - - distances = _get_boundary_point_distances(pixel_order, pixel, step) - margin_threshold_deg = margin_threshold / 3600.0 - - num_bounds = step * 4 - ra_chunks = np.array_split(r_asc, chunk_size) - dec_chunks = np.array_split(dec, chunk_size) - index = 0 - for chunk_index, ra_chunk in enumerate(ra_chunks): - dec_chunk = dec_chunks[chunk_index] - - chunk_len = len(ra_chunk) - - if chunk_len > 0: - ra_matrix = np.repeat(np.reshape([ra_chunk], (-1, 1)), num_bounds, axis=1) - dec_matrix = np.repeat(np.reshape([dec_chunk], (-1, 1)), num_bounds, axis=1) - - bounds_ra_matrix = np.repeat(np.array([bounds[0]]), len(ra_chunk), axis=0) - bounds_dec_matrix = np.repeat(np.array([bounds[1]]), len(dec_chunk), axis=0) - - points_coords = SkyCoord(ra=ra_matrix, dec=dec_matrix, unit="deg") - bounds_coords = SkyCoord(ra=bounds_ra_matrix, dec=bounds_dec_matrix, unit="deg") - - separations = points_coords.separation(bounds_coords).deg - - bisectors = np.apply_along_axis( - _find_minimum_distance, - 1, - separations, - distances=distances, - margin_threshold=margin_threshold_deg, - ) - margin_check[index : index + chunk_len] = bisectors <= margin_threshold_deg - index += chunk_len - - del ra_matrix, dec_matrix, bounds_ra_matrix, bounds_dec_matrix - del points_coords, bounds_coords, separations - - del distances - - return margin_check - - -# numba jit compiler doesn't count for coverage tests, so we'll set no cover. -@njit -def _find_minimum_distance(separations, distances, margin_threshold): # pragma: no cover - """Find the minimum distance between a given datapoint and a healpixel""" - minimum_index = np.argmin(separations) - - if separations[minimum_index] <= margin_threshold: - return separations[minimum_index] - - plus_one, minus_one = minimum_index + 1, minimum_index - 1 - num_seps = len(separations) - if minimum_index == 0: - minus_one = num_seps - 1 - if minimum_index == num_seps - 1: - plus_one = 0 - - if separations[minus_one] <= separations[plus_one]: - other_index = minus_one - else: - other_index = plus_one - - side_x = np.radians(separations[minimum_index]) - side_y = np.radians(separations[other_index]) - if (0 in (minimum_index, other_index)) and (num_seps - 1 in (minimum_index, other_index)): - side_z = distances[-1] - else: - side_z = distances[min(minimum_index, other_index)] - side_z *= np.pi / 180.0 - - ang = np.arccos((np.cos(side_y) - (np.cos(side_x) * np.cos(side_z))) / (np.sin(side_x) * np.sin(side_z))) - bisector = np.degrees(np.arcsin(np.sin(ang) * np.sin(side_x))) - - return bisector - - -def _get_boundary_point_distances(order, pixel, step): - """Get the distance of the segments between healpixel boundary points""" - boundary_points = hp.vec2dir(hp.boundaries(2**order, pixel, step=step, nest=True), lonlat=True) - - # shift forward all the coordinates by one - shift_ra = np.roll(boundary_points[0], 1) - shift_dec = np.roll(boundary_points[1], 1) - - coord_set_1 = SkyCoord(ra=boundary_points[0], dec=boundary_points[1], unit="deg") - coord_set_2 = SkyCoord(ra=shift_ra, dec=shift_dec, unit="deg") - - separations = coord_set_1.separation(coord_set_2).degree - return separations diff --git a/tests/hats/pixel_math/test_margin_bounding.py b/tests/hats/pixel_math/test_margin_bounding.py deleted file mode 100644 index 98945f2d..00000000 --- a/tests/hats/pixel_math/test_margin_bounding.py +++ /dev/null @@ -1,109 +0,0 @@ -"""Tests of margin bounding utility functions""" - -import numpy as np -import numpy.testing as npt -import pytest - -import hats.pixel_math as pm - - -def test_check_margin_bounds(): - """Make sure check_margin_bounds works properly""" - - ras = np.array([42.4704538, 56.25, 56.25, 50.61225197]) - decs = np.array([1.4593925, 9.6, 10.0, 14.4767556]) - - checks = pm.check_margin_bounds(ras, decs, 3, 4, 360.0) - - expected = np.array([False, True, False, True]) - - npt.assert_array_equal(checks, expected) - - -def test_check_margin_bounds_low_order(): - """Make sure check_margin_bounds works when using a low order healpixel""" - - ras = np.array( - [ - 142.4704538, - 45.09, - 45.2, - 37.31343956517391, - 42.649354753311535, - 32.62796809350278, - 39.89468227954832, - 27.718121934039974, - ] - ) - - decs = np.array( - [ - 1.4593925, - 0.0, - 0.0, - 6.566326903165274, - 2.005185097251452, - 10.597884275167646, - 4.465967883812584, - 14.959672304191956, - ] - ) - - checks = pm.check_margin_bounds(ras, decs, 0, 4, 360.0) - - expected = np.array([False, True, False, True, True, True, True, True]) - - npt.assert_array_equal(checks, expected) - - -def test_check_polar_margin_bounds_north(): - """Make sure check_polar_margin_bounds works at the north pole""" - order = 0 - pix = 1 - - ras = np.array([89, -179, -45]) - decs = np.array([89.9, 89.9, 85.0]) - - vals = pm.check_margin_bounds(ras, decs, order, pix, 360.0) - - expected = np.array([True, True, False]) - - npt.assert_array_equal(vals, expected) - - -def test_check_polar_margin_bounds_south(): - """Make sure check_polar_margin_bounds works at the south pole""" - order = 0 - pix = 9 - - ras = np.array([89, -179, -45]) - decs = np.array([-89.9, -89.9, -85.0]) - - vals = pm.check_margin_bounds(ras, decs, order, pix, 360.0) - - expected = np.array([True, True, False]) - - npt.assert_array_equal(vals, expected) - - -def test_check_margin_bounds_uneven(): - """Ensure check_margin_bounds fails when ra and dec arrays are unbalanced.""" - - ras = np.array([42.4704538, 56.25, 56.25, 50.61225197]) - decs = np.array([1.4593925, 9.6, 10.0]) - - with pytest.raises(ValueError, match="length of r_asc"): - pm.check_margin_bounds(ras, decs, 3, 4, 360.0) - - -def test_check_margin_bounds_empty(): - """Ensure check_margin_bounds works when passed empty coordinate arrays.""" - - ras = np.array([]) - decs = np.array([]) - - checks = pm.check_margin_bounds(ras, decs, 3, 4, 360.0) - - expected = np.array([]) - - npt.assert_array_equal(checks, expected) diff --git a/tests/hats/pixel_math/test_spatial_index.py b/tests/hats/pixel_math/test_spatial_index.py index 587482f0..ca9b04eb 100644 --- a/tests/hats/pixel_math/test_spatial_index.py +++ b/tests/hats/pixel_math/test_spatial_index.py @@ -132,11 +132,11 @@ def test_spatial_index_to_healpix_low_order(): def test_healpix_to_spatial_index_single(): orders = [3, 3, 4, 1] pixels = [0, 12, 1231, 11] - pixels_at_high_order = [p * (4 ** (SPATIAL_INDEX_ORDER - o)) for o, p in zip(orders, pixels)] - lon, lat = hp.pix2ang( - [2**SPATIAL_INDEX_ORDER] * len(orders), pixels_at_high_order, nest=True, lonlat=True - ) - actual_spatial_indices = compute_spatial_index(lon, lat) + + ra = [45.0, 45.0, 0.0, 225.0] + dec = [7.11477952e-08, 1.94712207e01, 1.44775123e01, 4.18103150e01] + + actual_spatial_indices = compute_spatial_index(ra, dec) test_spatial_indices = [healpix_to_spatial_index(o, p) for o, p in zip(orders, pixels)] assert np.all(test_spatial_indices == actual_spatial_indices) @@ -144,10 +144,9 @@ def test_healpix_to_spatial_index_single(): def test_healpix_to_spatial_index_array(): orders = [3, 3, 4, 1] pixels = [0, 12, 1231, 11] - pixels_at_high_order = [p * (4 ** (SPATIAL_INDEX_ORDER - o)) for o, p in zip(orders, pixels)] - lon, lat = hp.pix2ang( - [2**SPATIAL_INDEX_ORDER] * len(orders), pixels_at_high_order, nest=True, lonlat=True - ) - actual_spatial_indices = compute_spatial_index(lon, lat) + + ra = [45.0, 45.0, 0.0, 225.0] + dec = [7.11477952e-08, 1.94712207e01, 1.44775123e01, 4.18103150e01] + actual_spatial_indices = compute_spatial_index(ra, dec) test_spatial_indices = healpix_to_spatial_index(orders, pixels) assert np.all(test_spatial_indices == actual_spatial_indices) From 2dca8ec18b1d867abd29b1ca340fd09a19d30301 Mon Sep 17 00:00:00 2001 From: Melissa DeLucchi Date: Mon, 18 Nov 2024 11:08:24 -0500 Subject: [PATCH 02/19] Remove unused unseen. --- src/hats/pixel_math/healpix_shim.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/src/hats/pixel_math/healpix_shim.py b/src/hats/pixel_math/healpix_shim.py index 44aed38b..bee05e47 100644 --- a/src/hats/pixel_math/healpix_shim.py +++ b/src/hats/pixel_math/healpix_shim.py @@ -40,10 +40,6 @@ def ring2nest(*args, **kwargs): return hp.ring2nest(*args, **kwargs) -def unseen_pixel(): - return hp.pixelfunc.UNSEEN - - ## Query From f732d786205e976b488eef42cecc77333f9e50ba Mon Sep 17 00:00:00 2001 From: Melissa DeLucchi <113376043+delucchi-cmu@users.noreply.github.com> Date: Mon, 18 Nov 2024 12:25:46 -0500 Subject: [PATCH 03/19] Remove margin fine filtering. (#421) * Remove margin fine filtering. * Remove unused unseen. --- .github/workflows/pre-commit-ci.yml | 2 +- .github/workflows/testing-and-coverage.yml | 2 +- src/hats/pixel_math/__init__.py | 1 - src/hats/pixel_math/healpix_shim.py | 24 --- src/hats/pixel_math/margin_bounding.py | 139 ------------------ tests/hats/pixel_math/test_margin_bounding.py | 109 -------------- tests/hats/pixel_math/test_spatial_index.py | 19 ++- 7 files changed, 11 insertions(+), 285 deletions(-) delete mode 100644 src/hats/pixel_math/margin_bounding.py delete mode 100644 tests/hats/pixel_math/test_margin_bounding.py diff --git a/.github/workflows/pre-commit-ci.yml b/.github/workflows/pre-commit-ci.yml index f31b3ee6..642e7b9d 100644 --- a/.github/workflows/pre-commit-ci.yml +++ b/.github/workflows/pre-commit-ci.yml @@ -7,7 +7,7 @@ on: push: branches: [ main ] pull_request: - branches: [ main ] + branches: [ main, margin ] jobs: pre-commit-ci: diff --git a/.github/workflows/testing-and-coverage.yml b/.github/workflows/testing-and-coverage.yml index 3d3b867f..8438b944 100644 --- a/.github/workflows/testing-and-coverage.yml +++ b/.github/workflows/testing-and-coverage.yml @@ -7,7 +7,7 @@ on: push: branches: [ main ] pull_request: - branches: [ main ] + branches: [ main, margin ] jobs: build: diff --git a/src/hats/pixel_math/__init__.py b/src/hats/pixel_math/__init__.py index 132eaa84..327d83f9 100644 --- a/src/hats/pixel_math/__init__.py +++ b/src/hats/pixel_math/__init__.py @@ -2,7 +2,6 @@ from .healpix_pixel import HealpixPixel from .healpix_pixel_convertor import HealpixInputTypes, get_healpix_pixel -from .margin_bounding import check_margin_bounds from .partition_stats import empty_histogram, generate_alignment, generate_histogram from .pixel_margins import get_margin from .spatial_index import compute_spatial_index, spatial_index_to_healpix diff --git a/src/hats/pixel_math/healpix_shim.py b/src/hats/pixel_math/healpix_shim.py index 766bb5d7..bee05e47 100644 --- a/src/hats/pixel_math/healpix_shim.py +++ b/src/hats/pixel_math/healpix_shim.py @@ -20,18 +20,10 @@ def npix2nside(param): return hp.npix2nside(param) -def nside2npix(param): - return hp.nside2npix(param) - - def order2npix(param): return hp.order2npix(param) -def npix2order(param): - return hp.npix2order(param) - - def nside2resol(*args, **kwargs): return hp.nside2resol(*args, **kwargs) @@ -48,10 +40,6 @@ def ring2nest(*args, **kwargs): return hp.ring2nest(*args, **kwargs) -def unseen_pixel(): - return hp.pixelfunc.UNSEEN - - ## Query @@ -63,10 +51,6 @@ def query_polygon(*args, **kwargs): return hp.query_polygon(*args, **kwargs) -def boundaries(*args, **kwargs): - return hp.boundaries(*args, **kwargs) - - def get_all_neighbours(*args, **kwargs): return hp.get_all_neighbours(*args, **kwargs) @@ -78,14 +62,6 @@ def ang2vec(*args, **kwargs): return hp.ang2vec(*args, **kwargs) -def pix2ang(*args, **kwargs): - return hp.pix2ang(*args, **kwargs) - - -def vec2dir(*args, **kwargs): - return hp.vec2dir(*args, **kwargs) - - def pix2xyf(*args, **kwargs): return hp.pix2xyf(*args, **kwargs) diff --git a/src/hats/pixel_math/margin_bounding.py b/src/hats/pixel_math/margin_bounding.py deleted file mode 100644 index 8234f740..00000000 --- a/src/hats/pixel_math/margin_bounding.py +++ /dev/null @@ -1,139 +0,0 @@ -"""Utilities to build bounding boxes around healpixels that include a neighor margin.""" - -import numpy as np -from astropy.coordinates import SkyCoord -from numba import njit - -import hats.pixel_math.healpix_shim as hp - - -def check_margin_bounds(r_asc, dec, pixel_order, pixel, margin_threshold, step=100, chunk_size=1000): - # pylint: disable=too-many-locals - """Check a set of points in ra and dec space to see if they fall within the - `margin_threshold` of a given healpixel. - - Args: - r_asc (:obj:`np.array`): one dimmensional array representing the ra of a - given set of points. - dec (:obj:`np.array`): one dimmensional array representing the dec of a - given set of points. - pixel_order (int): the order of the healpixel to be checked against. - pixel (int): the healpixel to be checked against. - margin_threshold (double): the size of the margin cache in arcseconds. - step (int): the amount of samples we'll get from the healpixel boundary - to test the datapoints against. the higher the step the more - granular the point checking and therefore the more accurate, however - using a smaller step value might be helpful for super large datasets - to save compute time if accuracy is less important. - chunk_size (int): the size of the chunk of data points we process on a single - thread at any given time. We only process a certain amount of the data - at once as trying to check all data points at once can lead to - memory issues. - Returns: - :obj:`numpy.array` of boolean values corresponding to the ra and dec - coordinates checked against whether a given point is within any of the - provided polygons. - """ - num_points = len(r_asc) - if num_points != len(dec): - raise ValueError( - f"length of r_asc ({num_points}) is different \ - from length of dec ({len(dec)})" - ) - - if num_points == 0: - return np.array([]) - - bounds = hp.vec2dir(hp.boundaries(2**pixel_order, pixel, step=step, nest=True), lonlat=True) - margin_check = np.full((num_points,), False) - - distances = _get_boundary_point_distances(pixel_order, pixel, step) - margin_threshold_deg = margin_threshold / 3600.0 - - num_bounds = step * 4 - ra_chunks = np.array_split(r_asc, chunk_size) - dec_chunks = np.array_split(dec, chunk_size) - index = 0 - for chunk_index, ra_chunk in enumerate(ra_chunks): - dec_chunk = dec_chunks[chunk_index] - - chunk_len = len(ra_chunk) - - if chunk_len > 0: - ra_matrix = np.repeat(np.reshape([ra_chunk], (-1, 1)), num_bounds, axis=1) - dec_matrix = np.repeat(np.reshape([dec_chunk], (-1, 1)), num_bounds, axis=1) - - bounds_ra_matrix = np.repeat(np.array([bounds[0]]), len(ra_chunk), axis=0) - bounds_dec_matrix = np.repeat(np.array([bounds[1]]), len(dec_chunk), axis=0) - - points_coords = SkyCoord(ra=ra_matrix, dec=dec_matrix, unit="deg") - bounds_coords = SkyCoord(ra=bounds_ra_matrix, dec=bounds_dec_matrix, unit="deg") - - separations = points_coords.separation(bounds_coords).deg - - bisectors = np.apply_along_axis( - _find_minimum_distance, - 1, - separations, - distances=distances, - margin_threshold=margin_threshold_deg, - ) - margin_check[index : index + chunk_len] = bisectors <= margin_threshold_deg - index += chunk_len - - del ra_matrix, dec_matrix, bounds_ra_matrix, bounds_dec_matrix - del points_coords, bounds_coords, separations - - del distances - - return margin_check - - -# numba jit compiler doesn't count for coverage tests, so we'll set no cover. -@njit -def _find_minimum_distance(separations, distances, margin_threshold): # pragma: no cover - """Find the minimum distance between a given datapoint and a healpixel""" - minimum_index = np.argmin(separations) - - if separations[minimum_index] <= margin_threshold: - return separations[minimum_index] - - plus_one, minus_one = minimum_index + 1, minimum_index - 1 - num_seps = len(separations) - if minimum_index == 0: - minus_one = num_seps - 1 - if minimum_index == num_seps - 1: - plus_one = 0 - - if separations[minus_one] <= separations[plus_one]: - other_index = minus_one - else: - other_index = plus_one - - side_x = np.radians(separations[minimum_index]) - side_y = np.radians(separations[other_index]) - if (0 in (minimum_index, other_index)) and (num_seps - 1 in (minimum_index, other_index)): - side_z = distances[-1] - else: - side_z = distances[min(minimum_index, other_index)] - side_z *= np.pi / 180.0 - - ang = np.arccos((np.cos(side_y) - (np.cos(side_x) * np.cos(side_z))) / (np.sin(side_x) * np.sin(side_z))) - bisector = np.degrees(np.arcsin(np.sin(ang) * np.sin(side_x))) - - return bisector - - -def _get_boundary_point_distances(order, pixel, step): - """Get the distance of the segments between healpixel boundary points""" - boundary_points = hp.vec2dir(hp.boundaries(2**order, pixel, step=step, nest=True), lonlat=True) - - # shift forward all the coordinates by one - shift_ra = np.roll(boundary_points[0], 1) - shift_dec = np.roll(boundary_points[1], 1) - - coord_set_1 = SkyCoord(ra=boundary_points[0], dec=boundary_points[1], unit="deg") - coord_set_2 = SkyCoord(ra=shift_ra, dec=shift_dec, unit="deg") - - separations = coord_set_1.separation(coord_set_2).degree - return separations diff --git a/tests/hats/pixel_math/test_margin_bounding.py b/tests/hats/pixel_math/test_margin_bounding.py deleted file mode 100644 index 98945f2d..00000000 --- a/tests/hats/pixel_math/test_margin_bounding.py +++ /dev/null @@ -1,109 +0,0 @@ -"""Tests of margin bounding utility functions""" - -import numpy as np -import numpy.testing as npt -import pytest - -import hats.pixel_math as pm - - -def test_check_margin_bounds(): - """Make sure check_margin_bounds works properly""" - - ras = np.array([42.4704538, 56.25, 56.25, 50.61225197]) - decs = np.array([1.4593925, 9.6, 10.0, 14.4767556]) - - checks = pm.check_margin_bounds(ras, decs, 3, 4, 360.0) - - expected = np.array([False, True, False, True]) - - npt.assert_array_equal(checks, expected) - - -def test_check_margin_bounds_low_order(): - """Make sure check_margin_bounds works when using a low order healpixel""" - - ras = np.array( - [ - 142.4704538, - 45.09, - 45.2, - 37.31343956517391, - 42.649354753311535, - 32.62796809350278, - 39.89468227954832, - 27.718121934039974, - ] - ) - - decs = np.array( - [ - 1.4593925, - 0.0, - 0.0, - 6.566326903165274, - 2.005185097251452, - 10.597884275167646, - 4.465967883812584, - 14.959672304191956, - ] - ) - - checks = pm.check_margin_bounds(ras, decs, 0, 4, 360.0) - - expected = np.array([False, True, False, True, True, True, True, True]) - - npt.assert_array_equal(checks, expected) - - -def test_check_polar_margin_bounds_north(): - """Make sure check_polar_margin_bounds works at the north pole""" - order = 0 - pix = 1 - - ras = np.array([89, -179, -45]) - decs = np.array([89.9, 89.9, 85.0]) - - vals = pm.check_margin_bounds(ras, decs, order, pix, 360.0) - - expected = np.array([True, True, False]) - - npt.assert_array_equal(vals, expected) - - -def test_check_polar_margin_bounds_south(): - """Make sure check_polar_margin_bounds works at the south pole""" - order = 0 - pix = 9 - - ras = np.array([89, -179, -45]) - decs = np.array([-89.9, -89.9, -85.0]) - - vals = pm.check_margin_bounds(ras, decs, order, pix, 360.0) - - expected = np.array([True, True, False]) - - npt.assert_array_equal(vals, expected) - - -def test_check_margin_bounds_uneven(): - """Ensure check_margin_bounds fails when ra and dec arrays are unbalanced.""" - - ras = np.array([42.4704538, 56.25, 56.25, 50.61225197]) - decs = np.array([1.4593925, 9.6, 10.0]) - - with pytest.raises(ValueError, match="length of r_asc"): - pm.check_margin_bounds(ras, decs, 3, 4, 360.0) - - -def test_check_margin_bounds_empty(): - """Ensure check_margin_bounds works when passed empty coordinate arrays.""" - - ras = np.array([]) - decs = np.array([]) - - checks = pm.check_margin_bounds(ras, decs, 3, 4, 360.0) - - expected = np.array([]) - - npt.assert_array_equal(checks, expected) diff --git a/tests/hats/pixel_math/test_spatial_index.py b/tests/hats/pixel_math/test_spatial_index.py index 587482f0..ca9b04eb 100644 --- a/tests/hats/pixel_math/test_spatial_index.py +++ b/tests/hats/pixel_math/test_spatial_index.py @@ -132,11 +132,11 @@ def test_spatial_index_to_healpix_low_order(): def test_healpix_to_spatial_index_single(): orders = [3, 3, 4, 1] pixels = [0, 12, 1231, 11] - pixels_at_high_order = [p * (4 ** (SPATIAL_INDEX_ORDER - o)) for o, p in zip(orders, pixels)] - lon, lat = hp.pix2ang( - [2**SPATIAL_INDEX_ORDER] * len(orders), pixels_at_high_order, nest=True, lonlat=True - ) - actual_spatial_indices = compute_spatial_index(lon, lat) + + ra = [45.0, 45.0, 0.0, 225.0] + dec = [7.11477952e-08, 1.94712207e01, 1.44775123e01, 4.18103150e01] + + actual_spatial_indices = compute_spatial_index(ra, dec) test_spatial_indices = [healpix_to_spatial_index(o, p) for o, p in zip(orders, pixels)] assert np.all(test_spatial_indices == actual_spatial_indices) @@ -144,10 +144,9 @@ def test_healpix_to_spatial_index_single(): def test_healpix_to_spatial_index_array(): orders = [3, 3, 4, 1] pixels = [0, 12, 1231, 11] - pixels_at_high_order = [p * (4 ** (SPATIAL_INDEX_ORDER - o)) for o, p in zip(orders, pixels)] - lon, lat = hp.pix2ang( - [2**SPATIAL_INDEX_ORDER] * len(orders), pixels_at_high_order, nest=True, lonlat=True - ) - actual_spatial_indices = compute_spatial_index(lon, lat) + + ra = [45.0, 45.0, 0.0, 225.0] + dec = [7.11477952e-08, 1.94712207e01, 1.44775123e01, 4.18103150e01] + actual_spatial_indices = compute_spatial_index(ra, dec) test_spatial_indices = healpix_to_spatial_index(orders, pixels) assert np.all(test_spatial_indices == actual_spatial_indices) From 3606aa62590f5f503c0c8eb6ad2992862a60412d Mon Sep 17 00:00:00 2001 From: Melissa DeLucchi Date: Mon, 18 Nov 2024 13:12:33 -0500 Subject: [PATCH 04/19] Use cdshealpix for margin pixel finding. --- docs/guide/pixel_math.md | 194 -------------------- docs/index.rst | 1 - src/hats/loaders/read_hats.py | 2 +- src/hats/pixel_math/healpix_shim.py | 16 +- src/hats/pixel_math/pixel_margins.py | 116 ++---------- tests/hats/catalog/test_catalog.py | 2 +- tests/hats/pixel_math/test_pixel_margins.py | 35 +--- 7 files changed, 19 insertions(+), 347 deletions(-) delete mode 100644 docs/guide/pixel_math.md diff --git a/docs/guide/pixel_math.md b/docs/guide/pixel_math.md deleted file mode 100644 index d02f90fb..00000000 --- a/docs/guide/pixel_math.md +++ /dev/null @@ -1,194 +0,0 @@ -# Pixel Math Appendix -A reference document for the various utility functions of `hats/pixel_math`. - -## Pixel Margins -The functions made to find the pixels that make up the border region of a given -healpixel. Primarly used as a way to speed up the neighbor/margin caching code -for [hats-import](https://github.com/astronomy-commons/hats-import/). Code -originally created by Mario Juric for HIPS, found -[here](https://github.com/mjuric/HIPS/blob/feature/multiprocess/hipscat/healpix.py). - -### get_edge -Given a pixel pix at some order, return all -pixels order dk _higher_ than pix's order that line -pix's edge (or a corner). - -pix: the pixel ID for which the margin is requested -dk: the requested order of edge pixel IDs, relative to order of pix -edge: which edge/corner to return (NE edge=0, E corner=1, SE edge = 2, ....) - -#### Algorithm - -If you look at how the NEST indexing scheme works, a pixel at some order is -subdivided into four subpixel at every subsequent order in such a way that the south -subpixel index equals `4*pix`, east is `4*pix+1`, west is `4*pix+2`, north is `4*pix+3`: - -``` - 4*pix+3 -pix -> 4*pix+2 4*pix+1 - 4*pix -``` - -Further subdivisions split up each of those sub pixels accordingly. For example, -the eastern subpixel (`4*pix+1`) gets divided up into four more: - -``` -S=4*(4*pix+1), E=4*(4*pix+1)+1, W=4*(4*pix+1)+2, N=4*(4*pix+1)+3 - - 4*(4*pix+3)+3 - 4*pix+3 4*(4*pix+3)+2 4*(4*pix+3)+1 - 4*(4*pix+2)+3 4*(4*pix+3) 4*(4*pix+1)+3 -pix ===> 4*pix+2 4*pix+1 ===> 4*(4*pix+2)+2 4*(4*pix+2)+1 4*(4*pix+1)+2 4*(4*pix+1)+1 - 4*(4*pix+2) 4*(4*pix)+3 4*(4*pix+1) - 4*pix 4*(4*pix)+2 4*(4*pix)+1 - 4*(4*pix) -``` -etcetera... - -We can see that the edge indices follow a pattern. For example, after two -subdivisions the south-east edge would consist of pixels: -``` -4*(4*pix), 4*(4*pix)+1 -4*(4*pix+1), 4*(4*pix+1)+1 -``` -with the top coming from subdividing the southern, and bottom the eastern pixel. -and so on, recursively, for more subdivisions. Similar patterns are identifiable -with other edges. - -This can be compactly written as: - -``` -4*(4*(4*pix + i) + j) + k .... -``` - -with i, j, k, ... being indices whose choice specifies which edge we get. -For example iterating through, i = {0, 1}, j = {0, 1}, k = {0, 1} generates -indices for the 8 pixels of the south-east edge for 3 subdivisions. Similarly, -for the north-west edge the index values would loop through {2, 3}, etc. - -This can be expanded as: - -``` -4**dk*pix + 4**(dk-1)*i + 4**(dk-2)*j + 4**(dk-3) * k + ... -``` - -where dk is the number of subdivisions. E.g., for dk=3, this would equal: - -``` -4**3*pix + 4**2*i + 4**1*j + k -``` - -When written with bit-shift operators, another interpretation becomes clearer: - -``` -pix << 6 + i << 4 + j << 2 + k -``` - -or if we look at the binary representation of the resulting number: - -``` - [pix][ii][jj][kk] -``` - -Where ii, jj, kk are the bits of _each_ of the i,j,k indices. That is, to get the -list of subpixels on the edge, we bit-shift the base index pix by 2*dk to the left, -and fill in the rest with all possible combinations of ii, jj, kk indices. For example, -the northeast edge has index values {2, 3} = {0b10, 0b11}, so for (say) pix=8=0b1000, the -northeast edge indices after two subdivisions are equal to: - -``` -0b1000 10 10 = 138 -0b1000 10 11 = 139 -0b1000 11 10 = 148 -0b1000 11 11 = 143 -``` - -Note that for a given edge and dk, the suffix of each edge index does not depend on -the value of pix (pix is bit-shifted to the right and added). This means it can be -precomputed & stored for repeated reuse. - -#### Implementation - -The implementation is based on the binary digit interpretation above. For a requested -edge and dk subdivision level, we generate (and cache) the suffixes. Then, for a given -pixel pix, the edge pixels are readily computed by left-shifting pix by 2*dk places, -and summing (== or-ing) it with the suffix array. - -### get_margin -Returns the list of indices of pixels of order (kk+dk) that -border the pixel pix of order kk. - -#### Algorithm -Given a pixel pix, find all of its neighbors. Then find the -edge at level (kk+dk) for each of the neighbors. - -This is relatively straightforward in the equatorial faces of the Healpix -sphere -- e.g., one takes the SW neighbor and requests its NE edge to get -the margin on that side of the pixel, then the E corner of the W neighbor, -etc... - -This breaks down on pixels that are at the edge of polar faces. There, -the neighbor's sense of direction _rotates_ when switching from face to -face. For example at order=2, pixel 5 is bordered by pixel 26 to the -northeast (see the Fig3, bottom, https://healpix.jpl.nasa.gov/html/intronode4.htm). -But to pixel 5 that is the **northwest** edge (and not southwest, as you'd -expect; run `hp.get_all_neighbours(4, 26, nest=True)` and -`hp.get_all_neighbours(4, 5, nest=True)` to verify these claims). - -This is because directions rotate 90deg clockwise when moving from one -polar face to another in easterly direction (e.g., from face 0 to face 1). -We have to identify when this happens, and change the edges we request -for such neighbors. Mathematically, this rotation is equal to adding +2 -(modulo 8) to the requested edge in get_edge(). I.e., for the -pixel 5 example, rather than requesting SE=4 edge of pixel 26, -we request SE+2=6=NE edge (see the comments in the definition of _edge_vectors -near get_edge()). - -This index shift generalizes to `2 * (neighbor_face - pix_face)` for the case -when _both_ faces are around the pole (either north or south). In the south, -the rotation is in the opposite direction (ccw) -- so the sign of the shift -changes. The generalized formula also handles the pixels whose one vertex -is the pole (where all three neighbors sharing that vertex are on different -faces). - -#### Implementation -Pretty straightforward implementation of the algorithm above. - -#### Algorithm -In the ring id scheme for `healpix`, the first and last 4 pixels are the polar pixels. To determine whether a nest scheme pixel is at the poles, all we have to do is convert the pixel into the ring scheme and determine if it falls at the beginning or end of the range 0 -> npix. So, if in the ring scheme the pix is `<= 3`, or `>= npix - 4`, we can safely assume that it is a polar pixel. - -## Margin Bounding -After constraining our data points using the `get_margin` code in `pixel_margins`, we then move on to our more accurate bound checking by calculating the distance between a given data point and the approximate boundaries of the healpixel. - -![diagram showing the spherical triangle formed by the data point and two boundary points with a perpendicular bisector.](margin_cache_diagram.png) - -To approximate this distance, we get the perpendicular bisector of the triangle made with the datapoint and the two closest boundary points. To find the the perpendicular bisector of this triangle (which is approximately the shortest distance between the data point and the healpix polygon), we will have to find - -$\sin(a) = \sin(A)\sin(x)$ - -$a = \arcsin(\sin(A),\sin(x))$ - -where $a$ is the bisector, $x$ is the arc between our data point and one of the boundary points, and $A$ is the angle between that arc and the arc of the healpix boundary. - -This identity comes from [Napier's rules for right spherical triangles](https://en.wikipedia.org/wiki/Spherical_trigonometry#Napier's_rules_for_right_spherical_triangles). - -we already have $x$ from our separation calculations, but we're going to need to calculate the angle $A$. to do this, we can make use of the [spherical law of cosines](https://en.wikipedia.org/wiki/Spherical_law_of_cosines). let's use the identity - -$\cos(y) = \cos(x)\cos(z) + \sin(x)\sin(z)\cos(A)$ - -from there we can solve for A with - -$\sin(x)\sin(z)\cos(A) = \cos(y) - \cos(x)\cos(z)$ - -$\cos(A) = {\cos(y) - \cos(x)\cos(z) \over \sin(x)\sin(z)}$ - -$A = {\arccos({\cos(y) - \cos(x)\cos(z) \over \sin(x)\sin(z)})}$ - -and now we have all the information to find the bisector and compare it to our margin threshold! -### check_margin_bounds -Given a set of ra and dec coordinates as well as the healpixel, healpix order, margin threshold, and the step value of sampled healpix boundaries, return a 1-dimmensional array of booleans on whether a given ra and dec coordinate pair are contained within any of the given bounding boxes. - -#### Implementation -We take a chunk of our points and generate a set matrix of points with a shape (size_of_chunk, number_of_boundary_points) with all of our data points repeated along the first axis. Then, we generate a matrix of the same shape with the boundary points repeated as well and then put both of those matrices into a `SkyCoord` so that we can calculate the separations between all of them. We also find all of the distances between adjacent boundary points, so that we can find know the distance of each boundary segment. For each data point, we find the boundary point with the minimum separation, and then find which of the two adjacent boundary points has the smaller separation from the datapoint. From there, we can do the above calculations to get the perpendicular bisector and compare this against the `margin_threshold`. For checking around the corners of the healpixel, where the data point may not actually form a neat triangle with it's closest boundary points, we check to see if the minimum distance is less than the margin threshold and return True if it is. This way, we'll keep any data points that fall within the margin threshold but would return NaN with our perpendicular bisector calculations. - -To speed up our calculations, the inner loops of calculations is compiled with the `numba` JIT compiler. diff --git a/docs/index.rst b/docs/index.rst index e3228de9..8f43a115 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -36,7 +36,6 @@ operations on top of these utilities. Some known extensions: :caption: Developers guide/contributing - Pixel math API Reference Getting Started diff --git a/src/hats/loaders/read_hats.py b/src/hats/loaders/read_hats.py index f071e2d2..15062036 100644 --- a/src/hats/loaders/read_hats.py +++ b/src/hats/loaders/read_hats.py @@ -72,7 +72,7 @@ def _read_moc_from_point_map(catalog_base_dir: str | Path | UPath) -> MOC | None if not file_io.does_file_or_directory_exist(point_map_path): return None fits_image = file_io.read_fits_image(point_map_path) - order = hp.nside2order(hp.npix2nside(len(fits_image))) + order = hp.npix2order(len(fits_image)) boolean_skymap = fits_image.astype(bool) ipix = np.where(boolean_skymap)[0] orders = np.full(ipix.shape, order) diff --git a/src/hats/pixel_math/healpix_shim.py b/src/hats/pixel_math/healpix_shim.py index bee05e47..9e23858e 100644 --- a/src/hats/pixel_math/healpix_shim.py +++ b/src/hats/pixel_math/healpix_shim.py @@ -8,18 +8,14 @@ ## Arithmetic conversions -def nside2order(param): - return hp.nside2order(param) +def npix2order(param): + return hp.nside2order(hp.npix2nside(param)) def order2nside(param): return hp.order2nside(param) -def npix2nside(param): - return hp.npix2nside(param) - - def order2npix(param): return hp.order2npix(param) @@ -51,10 +47,6 @@ def query_polygon(*args, **kwargs): return hp.query_polygon(*args, **kwargs) -def get_all_neighbours(*args, **kwargs): - return hp.get_all_neighbours(*args, **kwargs) - - ## Coordinate conversion @@ -62,10 +54,6 @@ def ang2vec(*args, **kwargs): return hp.ang2vec(*args, **kwargs) -def pix2xyf(*args, **kwargs): - return hp.pix2xyf(*args, **kwargs) - - ## FITS def read_map(*args, **kwargs): return hp.read_map(*args, **kwargs) diff --git a/src/hats/pixel_math/pixel_margins.py b/src/hats/pixel_math/pixel_margins.py index ef4298a4..d6200960 100644 --- a/src/hats/pixel_math/pixel_margins.py +++ b/src/hats/pixel_math/pixel_margins.py @@ -1,116 +1,24 @@ """Utilities for find the pixels of higher orders that surround a given healpixel.""" -import numpy as np +import cdshealpix -import hats.pixel_math.healpix_shim as hp -# see the documentation for get_edge() about this variable -_edge_vectors = [ - np.asarray([1, 3]), # NE edge - np.asarray([1]), # E corner - np.asarray([0, 1]), # SE edge - np.asarray([0]), # S corner - np.asarray([0, 2]), # SW edge - np.asarray([2]), # W corner - np.asarray([2, 3]), # NW edge - np.asarray([3]), # N corner -] - -# cache used by get_margin() -_suffixes = {} - - -def get_edge(d_order, pix, edge): - """Get all the pixels at order kk+dk bordering pixel pix. - See hats/pixel_math/README.md for more info. +def get_margin(order, pixel, delta_order): + """Get all the pixels at order order+delta_order bordering pixel pixel. Args: - dk (int): the change in k that we wish to find the margins for. + order (int): the healpix order of pixel. pix (int): the healpix pixel to find margin pixels of. - edge (int): the edge we want to find for the given pixel. (0-7) - - - 0: NE edge - - 1: E corner - - 2: SE edge - - 3: S corner - - 4: SW edge - - 5: W corner - - 6: NW edge - - 7: N corner + delta_order (int): the change in order that we wish to find the margins for. Returns: - one-dimensional numpy array of long integers, filled with the healpix pixels - at order kk+dk that border the given edge of pix. + one-dimensional numpy array of integers, filled with the healpix pixels + at order `order+delta_order` that border pixel. """ - if edge < 0 or edge > 7: - raise ValueError("edge can only be values between 0 and 7 (see docstring)") - - if (edge, d_order) in _suffixes: - pixel_edge = _suffixes[(edge, d_order)] - else: - # generate and cache the suffix: - - # generate all combinations of i,j,k,... suffixes for the requested edge - # See https://stackoverflow.com/a/35608701 - grid = np.array(np.meshgrid(*[_edge_vectors[edge]] * d_order)) - pixel_edge = grid.T.reshape(-1, d_order) - # bit-shift each suffix by the required number of bits - pixel_edge <<= np.arange(2 * (d_order - 1), -2, -2) - # sum them up row-wise, generating the suffix - pixel_edge = pixel_edge.sum(axis=1) - # cache for further reuse - _suffixes[(edge, d_order)] = pixel_edge - - # append the 'pix' preffix - pixel_edge = (pix << 2 * d_order) + pixel_edge - return pixel_edge - - -def get_margin(order, pix, d_order): - """Get all the pixels at order order+dk bordering pixel pix. - See hats/pixel_math/README.md for more info. - - Args: - order (int): the healpix order of pix. - pix (int): the healpix pixel to find margin pixels of. - d_order (int): the change in k that we wish to find the margins for. - Must be greater than kk. - Returns: - one-dimensional numpy array of long integers, filled with the healpix pixels - at order kk+dk that border pix. - """ - if d_order < 1: - raise ValueError("dk must be greater than order") - - nside = hp.order2nside(order) - - # get all neighboring pixels - neighbors = hp.get_all_neighbours(nside, pix, nest=True) - - # get the healpix faces IDs of pix and the neighboring pixels - _, _, pix_face = hp.pix2xyf(nside, pix, nest=True) - _, _, faces = hp.pix2xyf(nside, neighbors, nest=True) - - # indices which tell get_edge() which edge/vertex to return - # for a given pixel. The default order is compatible with the - # order returned by hp.get_all_neighbours(). - which = np.arange(8) - - # northern hemisphere; 90deg cw rotation for every +1 face increment - if pix_face < 4: - mask = faces < 4 - which[mask] += 2 * (faces - pix_face)[mask] - which %= 8 - # southern hemisphere; 90deg ccw rotation for every +1 face increment - elif pix_face >= 8: - mask = faces >= 8 - which[mask] -= 2 * (faces - pix_face)[mask] - which %= 8 + edges, corners = cdshealpix.external_neighbours(ipix=pixel, depth=order, delta_depth=delta_order) - # get all edges/vertices - # (making sure we skip -1 entries, for pixels with seven neighbors) margins = [] - for edge, pixel in zip(which, neighbors): - if pixel != -1: - margins.append(get_edge(d_order, pixel, edge)) - margins = np.concatenate(margins) + margins.extend(edges[0]) + margins.extend(corners[0]) + margins.sort() + margins = [int(pix) for pix in margins if pix >= 0] return margins diff --git a/tests/hats/catalog/test_catalog.py b/tests/hats/catalog/test_catalog.py index 988e04a8..03883b2f 100644 --- a/tests/hats/catalog/test_catalog.py +++ b/tests/hats/catalog/test_catalog.py @@ -112,7 +112,7 @@ def test_load_catalog_small_sky_order1_moc(small_sky_order1_dir): assert isinstance(cat, Catalog) assert cat.moc is not None counts_skymap = read_fits_image(paths.get_point_map_file_pointer(small_sky_order1_dir)) - skymap_order = hp.nside2order(hp.npix2nside(len(counts_skymap))) + skymap_order = hp.npix2order(len(counts_skymap)) assert cat.moc.max_order == skymap_order assert np.all(cat.moc.flatten() == np.where(counts_skymap > 0)) diff --git a/tests/hats/pixel_math/test_pixel_margins.py b/tests/hats/pixel_math/test_pixel_margins.py index bea16ecb..5f7e8322 100644 --- a/tests/hats/pixel_math/test_pixel_margins.py +++ b/tests/hats/pixel_math/test_pixel_margins.py @@ -35,28 +35,13 @@ def test_get_margin(): 1119, ] ) + expected.sort() assert len(margins) == 20 npt.assert_array_equal(margins, expected) -def test_zero_dk(): - """Check that the code fails when trying to find margins at same order as the pixel.""" - with pytest.raises(ValueError) as value_error: - pm.get_margin(2, 2, 0) - - assert str(value_error.value) == "dk must be greater than order" - - -def test_negative_dk(): - """Check that the code fails when trying to find margins at a higher order than the pixel.""" - with pytest.raises(ValueError) as value_error: - pm.get_margin(2, 2, -1) - - assert str(value_error.value) == "dk must be greater than order" - - def test_polar_edge(): """Check that the code works when trying to find margins around the north pole.""" margins = pm.get_margin(2, 5, 2) @@ -84,6 +69,7 @@ def test_polar_edge(): 1519, ] ) + expected.sort() assert len(margins) == 19 @@ -118,22 +104,7 @@ def test_polar_edge_south(): 527, ] ) + expected.sort() assert len(margins) == 20 npt.assert_array_equal(margins, expected) - - -def test_edge_negative_value(): - """Check to make sure get_edge fails when passed a negative edge value.""" - with pytest.raises(ValueError) as value_error: - pm.get_edge(2, 5, -1) - - assert str(value_error.value) == "edge can only be values between 0 and 7 (see docstring)" - - -def test_edge_greater_than_7(): - """Check to make sure get_edge fails when passed an edge value greater than 7.""" - with pytest.raises(ValueError) as value_error: - pm.get_edge(2, 5, 8) - - assert str(value_error.value) == "edge can only be values between 0 and 7 (see docstring)" From 56d343c901884af2b624ccbd28e108962dc007b3 Mon Sep 17 00:00:00 2001 From: Melissa DeLucchi Date: Mon, 18 Nov 2024 13:25:08 -0500 Subject: [PATCH 05/19] Unused import. --- tests/hats/pixel_math/test_pixel_margins.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/hats/pixel_math/test_pixel_margins.py b/tests/hats/pixel_math/test_pixel_margins.py index 5f7e8322..47ff6759 100644 --- a/tests/hats/pixel_math/test_pixel_margins.py +++ b/tests/hats/pixel_math/test_pixel_margins.py @@ -2,7 +2,6 @@ import numpy as np import numpy.testing as npt -import pytest import hats.pixel_math.pixel_margins as pm From 82f98dd18f418a0161fa8ca065ca0921979cf5fc Mon Sep 17 00:00:00 2001 From: Sean McGuire Date: Fri, 15 Nov 2024 19:03:03 -0500 Subject: [PATCH 06/19] add color_by_order option and plot on existing axes --- src/hats/inspection/visualize_catalog.py | 118 +++++++++++------ .../hats/inspection/test_visualize_catalog.py | 125 ++++++++++++++++++ 2 files changed, 203 insertions(+), 40 deletions(-) diff --git a/src/hats/inspection/visualize_catalog.py b/src/hats/inspection/visualize_catalog.py index 06344c99..433d0a96 100644 --- a/src/hats/inspection/visualize_catalog.py +++ b/src/hats/inspection/visualize_catalog.py @@ -81,7 +81,9 @@ def plot_pixels(catalog: HealpixDataset, plot_title: str | None = None, **kwargs ) -def plot_pixel_list(pixels: List[HealpixPixel], plot_title: str = "", projection="MOL", **kwargs): +def plot_pixel_list( + pixels: List[HealpixPixel], plot_title: str = "", projection="MOL", color_by_order=True, **kwargs +): """Create a visual map of the pixel density of a list of pixels. Args: @@ -89,6 +91,7 @@ def plot_pixel_list(pixels: List[HealpixPixel], plot_title: str = "", projection plot_title (str): heading for the plot projection (str): The projection to use. Available projections listed at https://docs.astropy.org/en/stable/wcs/supported_projections.html + color_by_order (bool): Whether to color the pixels by their order. True by default. kwargs: Additional args to pass to `plot_healpix_map` """ orders = np.array([p.order for p in pixels]) @@ -98,13 +101,16 @@ def plot_pixel_list(pixels: List[HealpixPixel], plot_title: str = "", projection order_map, projection=projection, title=plot_title, ipix=ipix, depth=orders, cbar=False, **kwargs ) col = ax.collections[0] - col_array = col.get_array() - plt.colorbar( - col, - boundaries=np.arange(np.min(col_array) - 0.5, np.max(col_array) + 0.6, 1), - ticks=np.arange(np.min(col_array), np.max(col_array) + 1), - label="HEALPix Order", - ) + if color_by_order: + col_array = col.get_array() + plt.colorbar( + col, + boundaries=np.arange(np.min(col_array) - 0.5, np.max(col_array) + 0.6, 1), + ticks=np.arange(np.min(col_array), np.max(col_array) + 1), + label="HEALPix Order", + ) + else: + col.set(cmap=None, norm=None, array=None) return fig, ax @@ -175,26 +181,7 @@ def plot_moc( return fig, ax -def cull_to_fov(depth_ipix_d: Dict[int, Tuple[np.ndarray, np.ndarray]], wcs): - """Culls a mapping of ipix to values to pixels that are inside the plot window defined by a WCS - - Modified from mocpy.moc.plot.utils.build_plotting_moc - - Any pixels too small are merged to a lower order, with the map values within a lower order pixel being - sampled - - Args: - depth_ipix_d (Dict[int, Tuple[np.ndarray, np.ndarray]]): Map of HEALPix order to a tuple of 2 arrays - (the ipix array of pixel numbers in NESTED ordering, and the values of the pixels) - wcs (astropy.WCS): The wcs object with the plot's projection - - Returns: - A new map with the same datatype of depth_ipix_d, with any pixels outside the plot's FOV removed, - and any pixels too small merged with their map values subsampled. - """ - - depth_ipix_d = _merge_too_small_pixels(depth_ipix_d, wcs) - +def get_fov_moc_from_wcs(wcs: WCS) -> MOC | None: # Get the MOC delimiting the FOV polygon width_px = int(wcs.wcs.crpix[0] * 2.0) # Supposing the wcs is centered in the axis height_px = int(wcs.wcs.crpix[1] * 2.0) @@ -223,10 +210,37 @@ def cull_to_fov(depth_ipix_d: Dict[int, Tuple[np.ndarray, np.ndarray]], wcs): warnings.filterwarnings("default") if np.isnan(ra_deg).any() or np.isnan(dec_deg).any(): - return depth_ipix_d + return None # Create a rough MOC (depth=3 is sufficient) from the viewport moc_viewport = MOC.from_polygon_skycoord(viewport, max_depth=3) + return moc_viewport + + +def cull_to_fov(depth_ipix_d: Dict[int, Tuple[np.ndarray, np.ndarray]], wcs): + """Culls a mapping of ipix to values to pixels that are inside the plot window defined by a WCS + + Modified from mocpy.moc.plot.utils.build_plotting_moc + + Any pixels too small are merged to a lower order, with the map values within a lower order pixel being + sampled + + Args: + depth_ipix_d (Dict[int, Tuple[np.ndarray, np.ndarray]]): Map of HEALPix order to a tuple of 2 arrays + (the ipix array of pixel numbers in NESTED ordering, and the values of the pixels) + wcs (astropy.WCS): The wcs object with the plot's projection + + Returns: + A new map with the same datatype of depth_ipix_d, with any pixels outside the plot's FOV removed, + and any pixels too small merged with their map values subsampled. + """ + + depth_ipix_d = _merge_too_small_pixels(depth_ipix_d, wcs) + + moc_viewport = get_fov_moc_from_wcs(wcs) + + if moc_viewport is None: + return depth_ipix_d output_dict = {} @@ -510,18 +524,26 @@ def initialize_wcs_axes( frame_class (Type[BaseFrame] | None): The class of the frame for the WCSAxes to be initialized with. if the `ax` kwarg is used, this value is ignored (By Default uses EllipticalFrame for full sky projection. If FOV is set, RectangularFrame is used) - ax (WCSAxes | None): The matplotlib axes to plot onto. If None, an axes will be created to be used. If - specified, the axes must be an astropy WCSAxes, and the `wcs` parameter must be set with the WCS - object used in the axes. (Default: None) - fig (Figure | None): The matplotlib figure to add the axes to. If None, one will be created, unless - ax is specified (Default: None) + ax (WCSAxes | None): The matplotlib axes to plot onto. If None, the current axes will be used. + If the current axes is not the correct WCSAxes type, a new figure and axes will be created to be + used. If specified, the axes must be an astropy WCSAxes, and the `wcs` parameter will be ignored + and the wcs of the axes used. (Default: None) + fig (Figure | None): The matplotlib figure to add the axes to. If None, the current figure will be + used, unless ax is specified, in which case the figure of the ax will be used. If there is no + current figure, one will be created. (Default: None) kwargs: additional kwargs to pass to figure initialization """ if fig is None: if ax is not None: + # Use figure from axes if ax provided fig = ax.get_figure() else: - fig = plt.figure(**kwargs) + if len(plt.get_fignums()) == 0: + # create new figure if no existing figure + fig = plt.figure(**kwargs) + else: + # use current figure if exists + fig = plt.gcf() if frame_class is None and fov is None and wcs is None: frame_class = EllipticalFrame if fov is None: @@ -529,20 +551,36 @@ def initialize_wcs_axes( if center is None: center = SkyCoord(0, 0, unit="deg", frame="icrs") if ax is None: + if len(fig.axes) > 0: + ax = plt.gca() + if isinstance(ax, WCSAxes): + # Use current axes if axes exists in figure and current axes is correct type + wcs = ax.wcs + return fig, ax, wcs + else: + # Plot onto new axes on new figure if current axes is not correct type + warnings.warn( + "Current axes is not of correct type WCSAxes. A new figure and axes will be used." + ) + fig = plt.figure(**kwargs) if wcs is None: + # Initialize wcs with params if no WCS provided wcs = WCS( fig, fov=fov, center=center, coordsys="icrs", - rotation=Angle(0, u.degree), + rotation=Angle(0, u.deg), projection=projection, ).w + # If no axes provided and no valid current axes then create new axes with the right projection ax = fig.add_subplot(1, 1, 1, projection=wcs, frame_class=frame_class) - elif wcs is None: - raise ValueError( - "if ax is provided, wcs must also be provided with the projection used in initializing ax" - ) + elif not isinstance(ax, WCSAxes): + # Error if provided axes is of wrong type + raise ValueError("ax is not of correct type WCSAxes") + else: + # Use WCS from provided axes if provided axes is correct type + wcs = ax.wcs return fig, ax, wcs diff --git a/tests/hats/inspection/test_visualize_catalog.py b/tests/hats/inspection/test_visualize_catalog.py index 397c8162..95f94c8b 100644 --- a/tests/hats/inspection/test_visualize_catalog.py +++ b/tests/hats/inspection/test_visualize_catalog.py @@ -6,6 +6,7 @@ import pytest from astropy.coordinates import Angle, SkyCoord from astropy.visualization.wcsaxes.frame import EllipticalFrame, RectangularFrame +from matplotlib import colors from matplotlib.colors import LogNorm, Normalize from matplotlib.path import Path from matplotlib.pyplot import get_cmap @@ -715,6 +716,103 @@ def test_plot_kwargs(): np.testing.assert_array_equal(col.get_array(), pix_map) +def test_plot_existing_fig(): + order = 3 + length = 10 + ipix = np.arange(length) + pix_map = np.arange(length) + depth = np.full(length, fill_value=order) + fig = plt.figure() + assert len(fig.axes) == 0 + fig_ret, ax_ret = plot_healpix_map(pix_map, ipix=ipix, depth=depth) + assert fig is fig_ret + ax = fig.get_axes()[0] + assert ax is ax_ret + assert len(ax.collections) == 1 + col = ax.collections[0] + paths = col.get_paths() + assert len(paths) == length + wcs = WCS( + fig, + fov=DEFAULT_FOV, + center=DEFAULT_CENTER, + coordsys=DEFAULT_COORDSYS, + rotation=DEFAULT_ROTATION, + projection=DEFAULT_PROJECTION, + ).w + for path, ipix in zip(paths, np.arange(len(pix_map))): + verts, codes = compute_healpix_vertices(order, np.array([ipix]), wcs) + np.testing.assert_array_equal(path.vertices, verts) + np.testing.assert_array_equal(path.codes, codes) + np.testing.assert_array_equal(col.get_array(), pix_map) + + +def test_plot_existing_wcsaxes(): + order = 3 + length = 10 + ipix = np.arange(length) + pix_map = np.arange(length) + depth = np.full(length, fill_value=order) + fig = plt.figure() + wcs = WCS( + fig, + fov=DEFAULT_FOV, + center=DEFAULT_CENTER, + coordsys=DEFAULT_COORDSYS, + rotation=DEFAULT_ROTATION, + projection=DEFAULT_PROJECTION, + ).w + ax = fig.add_subplot(1, 1, 1, projection=wcs) + assert len(fig.axes) == 1 + assert len(ax.collections) == 0 + fig_ret, ax_ret = plot_healpix_map(pix_map, ipix=ipix, depth=depth) + assert fig is fig_ret + assert ax is ax_ret + assert len(ax.collections) == 1 + col = ax.collections[0] + paths = col.get_paths() + assert len(paths) == length + for path, ipix in zip(paths, np.arange(len(pix_map))): + verts, codes = compute_healpix_vertices(order, np.array([ipix]), wcs) + np.testing.assert_array_equal(path.vertices, verts) + np.testing.assert_array_equal(path.codes, codes) + np.testing.assert_array_equal(col.get_array(), pix_map) + + +def test_plot_existing_wrong_axes(): + order = 3 + length = 10 + ipix = np.arange(length) + pix_map = np.arange(length) + depth = np.full(length, fill_value=order) + fig = plt.figure() + ax = fig.add_subplot(1, 1, 1) + assert len(fig.axes) == 1 + assert len(ax.collections) == 0 + with pytest.warns(match="WCSAxes"): + fig_ret, ax_ret = plot_healpix_map(pix_map, ipix=ipix, depth=depth) + assert fig is not fig_ret + assert ax is not ax_ret + assert len(ax.collections) == 0 + assert len(ax_ret.collections) == 1 + col = ax_ret.collections[0] + paths = col.get_paths() + wcs = WCS( + fig_ret, + fov=DEFAULT_FOV, + center=DEFAULT_CENTER, + coordsys=DEFAULT_COORDSYS, + rotation=DEFAULT_ROTATION, + projection=DEFAULT_PROJECTION, + ).w + assert len(paths) == length + for path, ipix in zip(paths, np.arange(len(pix_map))): + verts, codes = compute_healpix_vertices(order, np.array([ipix]), wcs) + np.testing.assert_array_equal(path.vertices, verts) + np.testing.assert_array_equal(path.codes, codes) + np.testing.assert_array_equal(col.get_array(), pix_map) + + def test_catalog_plot(small_sky_order1_catalog): fig, ax = plot_pixels(small_sky_order1_catalog) pixels = sorted(small_sky_order1_catalog.get_healpix_pixels()) @@ -738,6 +836,33 @@ def test_catalog_plot(small_sky_order1_catalog): assert ax.get_title() == f"Catalog pixel density map - {small_sky_order1_catalog.catalog_name}" +def test_catalog_plot_no_color_by_order(small_sky_order1_catalog): + fc = "white" + ec = "black" + fig, ax = plot_pixels(small_sky_order1_catalog, color_by_order=False, facecolor=fc, edgecolor=ec) + pixels = sorted(small_sky_order1_catalog.get_healpix_pixels()) + col = ax.collections[0] + paths = col.get_paths() + assert len(paths) == len(pixels) + wcs = WCS( + fig, + fov=DEFAULT_FOV, + center=DEFAULT_CENTER, + coordsys=DEFAULT_COORDSYS, + rotation=DEFAULT_ROTATION, + projection=DEFAULT_PROJECTION, + ).w + for p, path in zip(pixels, paths): + step = 2 ** (3 - p.order) + verts, codes = compute_healpix_vertices(p.order, np.array([p.pixel]), wcs, step=step) + np.testing.assert_array_equal(path.vertices, verts) + np.testing.assert_array_equal(path.codes, codes) + assert col.get_array() is None + np.testing.assert_array_equal(col.get_facecolor()[0], colors.to_rgba(fc)) + np.testing.assert_array_equal(col.get_edgecolor()[0], colors.to_rgba(ec)) + assert ax.get_title() == f"Catalog pixel density map - {small_sky_order1_catalog.catalog_name}" + + def test_plot_moc(small_sky_order1_catalog): small_sky_order1_catalog.moc.fill = MagicMock() _, ax = plot_moc(small_sky_order1_catalog.moc) From 1b52544f2bb99eb933ef8cce0b3c42dc60d02f8c Mon Sep 17 00:00:00 2001 From: Sean McGuire Date: Mon, 18 Nov 2024 18:55:04 -0500 Subject: [PATCH 07/19] reset matplotlib after tests --- .../hats/inspection/test_visualize_catalog.py | 56 ++++++++++++------- 1 file changed, 37 insertions(+), 19 deletions(-) diff --git a/tests/hats/inspection/test_visualize_catalog.py b/tests/hats/inspection/test_visualize_catalog.py index 95f94c8b..8cfdf1a1 100644 --- a/tests/hats/inspection/test_visualize_catalog.py +++ b/tests/hats/inspection/test_visualize_catalog.py @@ -35,6 +35,12 @@ DEFAULT_PROJECTION = "MOL" +@pytest.fixture(autouse=True) +def reset_matplotlib(): + yield + plt.close("all") + + def test_healpix_vertices(): depth = 3 ipix = np.array([10, 11]) @@ -87,7 +93,7 @@ def test_plot_healpix_pixels(): pix_map = np.arange(length) depth = np.full(length, fill_value=order) fig, ax = plot_healpix_map(pix_map, ipix=ipix, depth=depth) - assert len(ax.collections) == 1 + assert len(ax.collections) > 0 col = ax.collections[0] paths = col.get_paths() assert len(paths) == length @@ -121,7 +127,7 @@ def test_plot_healpix_pixels_different_order(): pix_map = np.arange(length) depth = np.full(length, fill_value=order) fig, ax = plot_healpix_map(pix_map, ipix=ipix, depth=depth) - assert len(ax.collections) == 1 + assert len(ax.collections) > 0 col = ax.collections[0] paths = col.get_paths() assert len(paths) == length @@ -149,7 +155,7 @@ def test_order_0_pixel_plots_with_step(): pix_map = np.array([map_value]) depth = np.array([0]) fig, ax = plot_healpix_map(pix_map, ipix=ipix, depth=depth) - assert len(ax.collections) == 1 + assert len(ax.collections) > 0 col = ax.collections[0] paths = col.get_paths() length = 1 @@ -178,7 +184,7 @@ def test_edge_pixels_split_to_order_7(): pix_map = np.array([map_value]) depth = np.array([0]) fig, ax = plot_healpix_map(pix_map, ipix=ipix, depth=depth) - assert len(ax.collections) == 1 + assert len(ax.collections) > 0 # Generate a dictionary of pixel indices that have sides that align with the meridian at ra = 180deg, the # right edge of the plot @@ -355,7 +361,7 @@ def test_plot_healpix_map(): ipix = np.arange(12 * 4**order) pix_map = np.arange(12 * 4**order) fig, ax = plot_healpix_map(pix_map) - assert len(ax.collections) == 1 + assert len(ax.collections) > 0 col = ax.collections[0] paths = col.get_paths() wcs = WCS( @@ -402,7 +408,7 @@ def test_plot_wcs_params(): center=SkyCoord(10, 10, unit="deg", frame="icrs"), projection="AIT", ) - assert len(ax.collections) == 1 + assert len(ax.collections) > 0 col = ax.collections[0] paths = col.get_paths() assert len(paths) == length @@ -438,7 +444,7 @@ def test_plot_wcs_params_frame(): projection="AIT", frame_class=EllipticalFrame, ) - assert len(ax.collections) == 1 + assert len(ax.collections) > 0 col = ax.collections[0] paths = col.get_paths() assert len(paths) == length @@ -473,7 +479,7 @@ def test_plot_fov_culling(): center=SkyCoord(10, 10, unit="deg", frame="icrs"), projection="AIT", ) - assert len(ax.collections) == 1 + assert len(ax.collections) > 0 col = ax.collections[0] paths = col.get_paths() wcs = WCS( @@ -520,7 +526,7 @@ def test_plot_wcs(): ).w fig2, ax = plot_healpix_map(pix_map, ipix=ipix, depth=depth, fig=fig, wcs=wcs) assert fig2 is fig - assert len(ax.collections) == 1 + assert len(ax.collections) > 0 col = ax.collections[0] paths = col.get_paths() assert len(paths) == length @@ -552,7 +558,7 @@ def test_plot_wcs_and_ax(): fig2, ax2 = plot_healpix_map(pix_map, ipix=ipix, depth=depth, fig=fig, wcs=wcs, ax=ax) assert fig2 is fig assert ax2 is ax - assert len(ax.collections) == 1 + assert len(ax.collections) > 0 col = ax.collections[0] paths = col.get_paths() assert len(paths) == length @@ -580,8 +586,20 @@ def test_plot_ax_no_wcs(): projection="AIT", ).w ax = fig.add_subplot(1, 1, 1, projection=wcs, frame_class=EllipticalFrame) - with pytest.raises(ValueError): - plot_healpix_map(pix_map, ipix=ipix, depth=depth, fig=fig, ax=ax) + assert len(ax.collections) == 0 + fig2, ax2 = plot_healpix_map(pix_map, ipix=ipix, depth=depth, fig=fig, ax=ax) + assert fig2 is fig + assert ax2 is ax + assert len(ax.collections) > 0 + col = ax.collections[0] + paths = col.get_paths() + assert len(paths) == length + for path, ipix in zip(paths, np.arange(len(pix_map))): + verts, codes = compute_healpix_vertices(order, np.array([ipix]), wcs) + np.testing.assert_array_equal(path.vertices, verts) + np.testing.assert_array_equal(path.codes, codes) + np.testing.assert_array_equal(col.get_array(), pix_map) + assert ax.get_transform("icrs") is not None def test_plot_cmaps(): @@ -592,7 +610,7 @@ def test_plot_cmaps(): pix_map = np.arange(length) depth = np.full(length, fill_value=order) fig, ax = plot_healpix_map(pix_map, ipix=ipix, depth=depth, cmap=cmap_name) - assert len(ax.collections) == 1 + assert len(ax.collections) > 0 col = ax.collections[0] paths = col.get_paths() assert col.get_cmap() == get_cmap(cmap_name) @@ -617,7 +635,7 @@ def test_plot_cmaps(): pix_map = np.arange(length) depth = np.full(length, fill_value=order) fig, ax = plot_healpix_map(pix_map, ipix=ipix, depth=depth, cmap=get_cmap(cmap_name)) - assert len(ax.collections) == 1 + assert len(ax.collections) > 0 col = ax.collections[0] paths = col.get_paths() assert col.get_cmap() == get_cmap(cmap_name) @@ -639,7 +657,7 @@ def test_plot_norm(): pix_map = np.arange(length) depth = np.full(length, fill_value=order) fig, ax = plot_healpix_map(pix_map, ipix=ipix, depth=depth, norm=norm) - assert len(ax.collections) == 1 + assert len(ax.collections) > 0 col = ax.collections[0] paths = col.get_paths() assert col.norm == norm @@ -668,7 +686,7 @@ def test_plot_no_cbar(): pix_map = np.arange(length) depth = np.full(length, fill_value=order) fig, ax = plot_healpix_map(pix_map, ipix=ipix, depth=depth, cbar=False) - assert len(ax.collections) == 1 + assert len(ax.collections) > 0 col = ax.collections[0] assert col.colorbar is None paths = col.get_paths() @@ -696,7 +714,7 @@ def test_plot_kwargs(): pix_map = np.arange(length) depth = np.full(length, fill_value=order) fig, ax = plot_healpix_map(pix_map, ipix=ipix, depth=depth, label=label) - assert len(ax.collections) == 1 + assert len(ax.collections) > 0 col = ax.collections[0] assert col.get_label() == label paths = col.get_paths() @@ -728,7 +746,7 @@ def test_plot_existing_fig(): assert fig is fig_ret ax = fig.get_axes()[0] assert ax is ax_ret - assert len(ax.collections) == 1 + assert len(ax.collections) > 0 col = ax.collections[0] paths = col.get_paths() assert len(paths) == length @@ -768,7 +786,7 @@ def test_plot_existing_wcsaxes(): fig_ret, ax_ret = plot_healpix_map(pix_map, ipix=ipix, depth=depth) assert fig is fig_ret assert ax is ax_ret - assert len(ax.collections) == 1 + assert len(ax.collections) > 0 col = ax.collections[0] paths = col.get_paths() assert len(paths) == length From 097180291ca81802f8d998cd52317557d3641740 Mon Sep 17 00:00:00 2001 From: Sean McGuire Date: Mon, 18 Nov 2024 19:13:57 -0500 Subject: [PATCH 08/19] pylint --- src/hats/inspection/visualize_catalog.py | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/src/hats/inspection/visualize_catalog.py b/src/hats/inspection/visualize_catalog.py index 433d0a96..cce75e75 100644 --- a/src/hats/inspection/visualize_catalog.py +++ b/src/hats/inspection/visualize_catalog.py @@ -182,6 +182,16 @@ def plot_moc( def get_fov_moc_from_wcs(wcs: WCS) -> MOC | None: + """Returns a MOC that matches the plot window defined by a WCS + + Modified from mocpy.moc.plot.utils.build_plotting_moc + + Args: + wcs (astropy.WCS): The wcs object with the plot's projection + + Returns: + The moc which defines the area of the sky that would be visible in a WCSAxes with the given WCS + """ # Get the MOC delimiting the FOV polygon width_px = int(wcs.wcs.crpix[0] * 2.0) # Supposing the wcs is centered in the axis height_px = int(wcs.wcs.crpix[1] * 2.0) @@ -557,12 +567,9 @@ def initialize_wcs_axes( # Use current axes if axes exists in figure and current axes is correct type wcs = ax.wcs return fig, ax, wcs - else: - # Plot onto new axes on new figure if current axes is not correct type - warnings.warn( - "Current axes is not of correct type WCSAxes. A new figure and axes will be used." - ) - fig = plt.figure(**kwargs) + # Plot onto new axes on new figure if current axes is not correct type + warnings.warn("Current axes is not of correct type WCSAxes. A new figure and axes will be used.") + fig = plt.figure(**kwargs) if wcs is None: # Initialize wcs with params if no WCS provided wcs = WCS( From 8f8bb6bb10a965874028b27ec456f2630482ae0f Mon Sep 17 00:00:00 2001 From: Sandro Campos Date: Tue, 19 Nov 2024 12:56:24 -0500 Subject: [PATCH 09/19] Migrate ang2vec (#426) * Use astropy for ang2vec conversion * Allow for kwargs in ang2vec --- src/hats/pixel_math/box_filter.py | 2 +- src/hats/pixel_math/healpix_shim.py | 8 ++++++-- src/hats/pixel_math/validators.py | 2 +- tests/hats/pixel_math/test_healpix_shim.py | 13 +++++++++++++ 4 files changed, 21 insertions(+), 4 deletions(-) diff --git a/src/hats/pixel_math/box_filter.py b/src/hats/pixel_math/box_filter.py index ae33ee45..d6155709 100644 --- a/src/hats/pixel_math/box_filter.py +++ b/src/hats/pixel_math/box_filter.py @@ -111,7 +111,7 @@ def _get_pixels_for_subpolygons(polygons: List[List[tuple[float, float]]], order nside = hp.order2nside(order) all_polygon_pixels = [] for vertices in polygons: - vertices = hp.ang2vec(*np.array(vertices).T, lonlat=True) + 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)) diff --git a/src/hats/pixel_math/healpix_shim.py b/src/hats/pixel_math/healpix_shim.py index 9e23858e..1afccda2 100644 --- a/src/hats/pixel_math/healpix_shim.py +++ b/src/hats/pixel_math/healpix_shim.py @@ -1,7 +1,9 @@ from __future__ import annotations +import astropy.units as u import healpy as hp import numpy as np +from astropy.coordinates import SkyCoord # pylint: disable=missing-function-docstring @@ -50,8 +52,10 @@ def query_polygon(*args, **kwargs): ## Coordinate conversion -def ang2vec(*args, **kwargs): - return hp.ang2vec(*args, **kwargs) +def ang2vec(ra, dec, **kwargs) -> np.ndarray: + """Converts ra and dec to cartesian coordinates on the unit sphere""" + coords = SkyCoord(ra=ra * u.deg, dec=dec * u.deg, **kwargs).cartesian + return np.array([coords.x.value, coords.y.value, coords.z.value]).T ## FITS diff --git a/src/hats/pixel_math/validators.py b/src/hats/pixel_math/validators.py index 6e60e432..9f87d6f3 100644 --- a/src/hats/pixel_math/validators.py +++ b/src/hats/pixel_math/validators.py @@ -84,7 +84,7 @@ def check_polygon_is_valid(vertices: np.ndarray): Returns: True if polygon is valid, False otherwise. """ - vertices_xyz = hp.ang2vec(*vertices.T, lonlat=True) + vertices_xyz = hp.ang2vec(*vertices.T) n_vertices = len(vertices_xyz) flip = 0 for i in range(n_vertices): diff --git a/tests/hats/pixel_math/test_healpix_shim.py b/tests/hats/pixel_math/test_healpix_shim.py index b83c978f..00ddc1bd 100644 --- a/tests/hats/pixel_math/test_healpix_shim.py +++ b/tests/hats/pixel_math/test_healpix_shim.py @@ -37,3 +37,16 @@ def test_order2mindist(): assert_allclose(hps.order2mindist(orders), min_distances, rtol=1e-2) assert_allclose(hps.order2mindist(17), 0.01677, rtol=1e-2) + + +def test_ang2vec(): + """Tests conversion of ra/dec to unit cartesian vectors""" + ra = [230.14467816, 110.40507118, 9.41764689, 245.5553117] + dec = [38.78080888, 17.09584081, -28.6352765, 5.41803306] + ra_rad = np.radians(ra) + dec_rad = np.radians(dec) + x = np.cos(ra_rad) * np.cos(dec_rad) + y = np.sin(ra_rad) * np.cos(dec_rad) + z = np.sin(dec_rad) + actual = np.asarray([x, y, z]).T + assert_array_equal(actual, hps.ang2vec(ra, dec)) From 7c8dfb7245fb7a6514e326e295530bec3fba7ccd Mon Sep 17 00:00:00 2001 From: Sean McGuire Date: Tue, 19 Nov 2024 13:59:49 -0500 Subject: [PATCH 10/19] use mindist instead of resol for margin filter --- src/hats/catalog/margin_cache/margin_catalog.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/hats/catalog/margin_cache/margin_catalog.py b/src/hats/catalog/margin_cache/margin_catalog.py index 65858e47..da40d160 100644 --- a/src/hats/catalog/margin_cache/margin_catalog.py +++ b/src/hats/catalog/margin_cache/margin_catalog.py @@ -33,8 +33,8 @@ def filter_by_moc(self, moc: MOC) -> Self: pixel sizes. """ max_order = moc.max_order - max_order_size = hp.nside2resol(2**max_order, arcmin=True) - if self.catalog_info.margin_threshold > max_order_size * 60: + max_order_size_arcsec = hp.order2mindist(max_order) * 60 + if self.catalog_info.margin_threshold > max_order_size_arcsec: raise ValueError( f"Cannot Filter Margin: Margin size {self.catalog_info.margin_threshold} is " f"greater than the size of a pixel at the highest order {max_order}." From 45512fc7473908a1a3f3ebe6c5b8159be1c9ab4c Mon Sep 17 00:00:00 2001 From: Sean McGuire Date: Wed, 20 Nov 2024 18:38:55 -0500 Subject: [PATCH 11/19] Replace healpy in healpix_shim pixel math operations --- src/hats/pixel_math/healpix_shim.py | 74 ++++++++--- tests/hats/pixel_math/test_healpix_shim.py | 142 +++++++++++++++++++++ 2 files changed, 201 insertions(+), 15 deletions(-) diff --git a/src/hats/pixel_math/healpix_shim.py b/src/hats/pixel_math/healpix_shim.py index 1afccda2..1c8a06f9 100644 --- a/src/hats/pixel_math/healpix_shim.py +++ b/src/hats/pixel_math/healpix_shim.py @@ -1,41 +1,85 @@ from __future__ import annotations +import math + import astropy.units as u +import cdshealpix import healpy as hp import numpy as np -from astropy.coordinates import SkyCoord +from astropy.coordinates import SkyCoord, Longitude, Latitude # pylint: disable=missing-function-docstring ## Arithmetic conversions -def npix2order(param): - return hp.nside2order(hp.npix2nside(param)) +MAX_HEALPIX_ORDER = 29 + + +def is_order_valid(order: int) -> bool: + return 0 <= order <= MAX_HEALPIX_ORDER + + +def npix2order(npix: int) -> int: + if npix <= 0: + raise ValueError("Invalid value for npix") + order = int(math.log2(npix / 12)) >> 1 + if not is_order_valid(order) or not 12 * (1 << (2 * order)) == npix: + raise ValueError("Invalid value for npix") + return order + + +def order2nside(order: int) -> int: + if not is_order_valid(order): + raise ValueError("Invalid value for order") + return 1 << order + + +def order2npix(order: int) -> int: + if not is_order_valid(order): + raise ValueError("Invalid value for order") + return 12 * (1 << (2 * order)) + +def nside2resol(nside: int, arcmin: bool = False) -> float: + resol_rad = np.sqrt(nside2pixarea(nside)) -def order2nside(param): - return hp.order2nside(param) + if arcmin: + return np.rad2deg(resol_rad) * 60 + return resol_rad -def order2npix(param): - return hp.order2npix(param) +def nside2pixarea(nside: int, degrees: bool = False) -> float: + npix = 12 * nside * nside + pix_area_rad = 4 * np.pi / npix + if degrees: + return pix_area_rad * (180 / np.pi) * (180 / np.pi) + return pix_area_rad -def nside2resol(*args, **kwargs): - return hp.nside2resol(*args, **kwargs) +def ang2pix(nside: int, theta: float, phi: float, nest: bool = False, lonlat: bool = False) -> int: + if not nest: + raise NotImplementedError("RING order ang2pix not supported") + order = nside2order(nside) + if lonlat: + ra = Longitude(theta, unit="deg") + dec = Latitude(phi, unit="deg") + else: + ra = Longitude(phi, unit="rad") + dec = Latitude(np.pi / 2 - theta, unit="rad") -def nside2pixarea(*args, **kwargs): - return hp.nside2pixarea(*args, **kwargs) + return cdshealpix.lonlat_to_healpix(ra, dec, order) -def ang2pix(*args, **kwargs): - return hp.ang2pix(*args, **kwargs) +def nside2order(nside: int) -> int: + npix = 12 * nside * nside + return npix2order(npix) -def ring2nest(*args, **kwargs): - return hp.ring2nest(*args, **kwargs) +def ring2nest(nside: int, ipix: int) -> int: + order = nside2order(nside) + return cdshealpix.from_ring(ipix, order) ## Query diff --git a/tests/hats/pixel_math/test_healpix_shim.py b/tests/hats/pixel_math/test_healpix_shim.py index 00ddc1bd..bc12ca9e 100644 --- a/tests/hats/pixel_math/test_healpix_shim.py +++ b/tests/hats/pixel_math/test_healpix_shim.py @@ -1,4 +1,7 @@ +import cdshealpix import numpy as np +import pytest +from astropy.coordinates import Longitude, Latitude from numpy.testing import assert_allclose, assert_array_equal from hats.pixel_math import healpix_shim as hps @@ -50,3 +53,142 @@ def test_ang2vec(): z = np.sin(dec_rad) actual = np.asarray([x, y, z]).T assert_array_equal(actual, hps.ang2vec(ra, dec)) + + +def test_npix2order(): + orders = [0, 1, 5, 10, 20, 29] + npix = [12 * (4**order) for order in orders] + test_orders = [hps.npix2order(x) for x in npix] + assert test_orders == orders + + +def test_npix2order_invalid(): + npixs = [-10, 0, 11, 13, 47, 49, 100000, 100000000000000000] + for npix in npixs: + with pytest.raises(ValueError, match="Invalid"): + hps.npix2order(npix) + + +def test_order2nside(): + orders = [0, 1, 5, 10, 20, 29] + expected_nsides = [2**x for x in orders] + test_nsides = [hps.order2nside(o) for o in orders] + assert test_nsides == expected_nsides + + +def test_order2nside_invalid(): + orders = [-1000, -1, 30, 4000] + for order in orders: + with pytest.raises(ValueError, match="Invalid"): + hps.order2nside(order) + + +def test_order2npix(): + orders = [0, 1, 5, 10, 20, 29] + npix = [12 * (4**order) for order in orders] + test_npix = [hps.order2npix(o) for o in orders] + assert test_npix == npix + + +def test_order2npix_invalid(): + orders = [-1000, -1, 30, 4000] + for order in orders: + with pytest.raises(ValueError, match="Invalid"): + hps.order2npix(order) + + +def test_nside2pixarea(): + orders = [0, 1, 5, 10, 20, 29] + nsides = [2**x for x in orders] + npix = [12 * (4**order) for order in orders] + pix_area_expected = [4 * np.pi / x for x in npix] + pix_area_test = [hps.nside2pixarea(nside) for nside in nsides] + assert pix_area_test == pix_area_expected + + +def test_nside2pixarea_degrees(): + orders = [0, 1, 5, 10, 20, 29] + nsides = [2**x for x in orders] + npix = [12 * (4**order) for order in orders] + pix_area_expected = [np.rad2deg(np.rad2deg(4 * np.pi / x)) for x in npix] + pix_area_test = [hps.nside2pixarea(nside, degrees=True) for nside in nsides] + assert pix_area_test == pix_area_expected + + +def test_nside2resol(): + orders = [0, 1, 5, 10, 20, 29] + nsides = [2**x for x in orders] + resol_expected = [np.sqrt(hps.nside2pixarea(nside)) for nside in nsides] + resol_test = [hps.nside2resol(nside) for nside in nsides] + assert resol_test == resol_expected + + +def test_nside2resol_arcmin(): + orders = [0, 1, 5, 10, 20, 29] + nsides = [2**x for x in orders] + resol_expected = [np.rad2deg(np.sqrt(hps.nside2pixarea(nside))) * 60 for nside in nsides] + resol_test = [hps.nside2resol(nside, arcmin=True) for nside in nsides] + assert resol_test == resol_expected + + +def test_ang2pix(): + orders = [0, 1, 5, 10, 20, 29] + ras = np.arange(-180.0, 180.0, 10.0) + decs = np.arange(-90.0, 90.0, 180 // len(ras)) + colats = 90 - decs + for order in orders: + expected_pixels = cdshealpix.lonlat_to_healpix( + Longitude(ras, unit="deg"), Latitude(decs, unit="deg"), order + ) + pixels = hps.ang2pix(hps.order2nside(order), np.deg2rad(colats), np.deg2rad(ras), nest=True) + assert np.all(pixels == expected_pixels) + + +def test_ang2pix_lonlat(): + orders = [0, 1, 5, 10, 20, 29] + ras = np.arange(-180.0, 180.0, 10.0) + decs = np.arange(-90.0, 90.0, 180 // len(ras)) + for order in orders: + expected_pixels = cdshealpix.lonlat_to_healpix( + Longitude(ras, unit="deg"), Latitude(decs, unit="deg"), order + ) + pixels = hps.ang2pix(hps.order2nside(order), ras, decs, nest=True, lonlat=True) + assert np.all(pixels == expected_pixels) + + +def test_ang2pix_ring(): + order = 4 + ras = np.arange(-180.0, 180.0, 10.0) + decs = np.arange(-90.0, 90.0, 180 // len(ras)) + colats = 90 - decs + with pytest.raises(NotImplementedError, match="RING"): + hps.ang2pix(hps.order2nside(order), ras, decs, lonlat=True) + with pytest.raises(NotImplementedError, match="RING"): + hps.ang2pix(hps.order2nside(order), colats, ras) + + +def test_ang2pix_invalid(): + orders = [0, 1, 5, 10, 20, 29] + invalid_orders = [-1000, -1, 30, 40] + ras = np.arange(-4000.0, 1000.0, 100.0) + decs = np.arange(-1000.0, 1000.0, 2000.0 // len(ras)) + colats = 90 - decs + for order in invalid_orders: + with pytest.raises(ValueError, match="Invalid"): + hps.ang2pix(int(2**order), np.deg2rad(colats), np.deg2rad(ras), nest=True) + with pytest.raises(ValueError, match="Invalid"): + hps.ang2pix(int(2**order), ras, decs, nest=True, lonlat=True) + for order in orders: + with pytest.raises(ValueError, match="angle"): + hps.ang2pix(hps.order2nside(order), np.deg2rad(colats), np.deg2rad(ras), nest=True) + with pytest.raises(ValueError, match="angle"): + hps.ang2pix(hps.order2nside(order), ras, decs, nest=True, lonlat=True) + + +def test_ring2nest(): + orders = [0, 1, 5, 10, 20, 29] + for order in orders: + ipix = np.arange(0, 12 * (4**order), 12 * (4**order) // 10) + expected_ring_ipix = cdshealpix.from_ring(ipix, order) + test_ring_ipix = hps.ring2nest(2**order, ipix) + assert np.all(test_ring_ipix == expected_ring_ipix) From 312f55772a40c332fb5fdd508a77b7277057c386 Mon Sep 17 00:00:00 2001 From: Sean McGuire Date: Wed, 20 Nov 2024 18:40:26 -0500 Subject: [PATCH 12/19] isort --- src/hats/pixel_math/healpix_shim.py | 2 +- tests/hats/pixel_math/test_healpix_shim.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/hats/pixel_math/healpix_shim.py b/src/hats/pixel_math/healpix_shim.py index 1c8a06f9..f46c57bc 100644 --- a/src/hats/pixel_math/healpix_shim.py +++ b/src/hats/pixel_math/healpix_shim.py @@ -6,7 +6,7 @@ import cdshealpix import healpy as hp import numpy as np -from astropy.coordinates import SkyCoord, Longitude, Latitude +from astropy.coordinates import Latitude, Longitude, SkyCoord # pylint: disable=missing-function-docstring diff --git a/tests/hats/pixel_math/test_healpix_shim.py b/tests/hats/pixel_math/test_healpix_shim.py index bc12ca9e..dd13654a 100644 --- a/tests/hats/pixel_math/test_healpix_shim.py +++ b/tests/hats/pixel_math/test_healpix_shim.py @@ -1,7 +1,7 @@ import cdshealpix import numpy as np import pytest -from astropy.coordinates import Longitude, Latitude +from astropy.coordinates import Latitude, Longitude from numpy.testing import assert_allclose, assert_array_equal from hats.pixel_math import healpix_shim as hps From 420a0443f1d1fb39158baaafe378a25a366d1421 Mon Sep 17 00:00:00 2001 From: Sean McGuire Date: Wed, 20 Nov 2024 18:48:41 -0500 Subject: [PATCH 13/19] support np arrays --- src/hats/pixel_math/healpix_shim.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/hats/pixel_math/healpix_shim.py b/src/hats/pixel_math/healpix_shim.py index f46c57bc..965f83d1 100644 --- a/src/hats/pixel_math/healpix_shim.py +++ b/src/hats/pixel_math/healpix_shim.py @@ -17,7 +17,7 @@ def is_order_valid(order: int) -> bool: - return 0 <= order <= MAX_HEALPIX_ORDER + return np.all(0 <= order) and np.all(order <= MAX_HEALPIX_ORDER) def npix2order(npix: int) -> int: From c7738f06c67d6a03f5ff4031f135654296c9fd5c Mon Sep 17 00:00:00 2001 From: Sean McGuire Date: Fri, 22 Nov 2024 11:27:13 -0500 Subject: [PATCH 14/19] use order instead of nside --- src/hats/catalog/partition_info.py | 2 +- src/hats/pixel_math/box_filter.py | 2 +- src/hats/pixel_math/healpix_shim.py | 35 ++++------- src/hats/pixel_math/partition_stats.py | 6 +- src/hats/pixel_math/spatial_index.py | 2 +- tests/hats/pixel_math/test_healpix_shim.py | 67 ++++++--------------- tests/hats/pixel_math/test_spatial_index.py | 18 ++---- 7 files changed, 40 insertions(+), 92 deletions(-) diff --git a/src/hats/catalog/partition_info.py b/src/hats/catalog/partition_info.py index d073702d..5188cca0 100644 --- a/src/hats/catalog/partition_info.py +++ b/src/hats/catalog/partition_info.py @@ -274,7 +274,7 @@ def calculate_fractional_coverage(self): """Calculate what fraction of the sky is covered by partition tiles.""" pixel_orders = [p.order for p in self.pixel_list] cov_order, cov_count = np.unique(pixel_orders, return_counts=True) - area_by_order = [hp.nside2pixarea(hp.order2nside(order), degrees=True) for order in cov_order] + area_by_order = [hp.order2pixarea(order, degrees=True) for order in cov_order] # 41253 is the number of square degrees in a sphere # https://en.wikipedia.org/wiki/Square_degree return (area_by_order * cov_count).sum() / (360**2 / np.pi) diff --git a/src/hats/pixel_math/box_filter.py b/src/hats/pixel_math/box_filter.py index d6155709..58c4b868 100644 --- a/src/hats/pixel_math/box_filter.py +++ b/src/hats/pixel_math/box_filter.py @@ -87,7 +87,7 @@ def _generate_dec_strip_moc(dec_range: Tuple[float, float], order: int) -> MOC: 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) + 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) diff --git a/src/hats/pixel_math/healpix_shim.py b/src/hats/pixel_math/healpix_shim.py index 965f83d1..96eb18b7 100644 --- a/src/hats/pixel_math/healpix_shim.py +++ b/src/hats/pixel_math/healpix_shim.py @@ -41,8 +41,8 @@ def order2npix(order: int) -> int: return 12 * (1 << (2 * order)) -def nside2resol(nside: int, arcmin: bool = False) -> float: - resol_rad = np.sqrt(nside2pixarea(nside)) +def order2resol(order: int, arcmin: bool = False) -> float: + resol_rad = np.sqrt(order2pixarea(order)) if arcmin: return np.rad2deg(resol_rad) * 60 @@ -50,35 +50,25 @@ def nside2resol(nside: int, arcmin: bool = False) -> float: return resol_rad -def nside2pixarea(nside: int, degrees: bool = False) -> float: - npix = 12 * nside * nside +def order2pixarea(order: int, degrees: bool = False) -> float: + npix = order2npix(order) pix_area_rad = 4 * np.pi / npix if degrees: return pix_area_rad * (180 / np.pi) * (180 / np.pi) return pix_area_rad -def ang2pix(nside: int, theta: float, phi: float, nest: bool = False, lonlat: bool = False) -> int: - if not nest: - raise NotImplementedError("RING order ang2pix not supported") - order = nside2order(nside) - if lonlat: - ra = Longitude(theta, unit="deg") - dec = Latitude(phi, unit="deg") - else: - ra = Longitude(phi, unit="rad") - dec = Latitude(np.pi / 2 - theta, unit="rad") - - return cdshealpix.lonlat_to_healpix(ra, dec, order) +def radec2pix(order: int, ra: float, dec: float) -> int: + if not is_order_valid(order): + raise ValueError("Invalid value for order") + ra = Longitude(ra, unit="deg") + dec = Latitude(dec, unit="deg") -def nside2order(nside: int) -> int: - npix = 12 * nside * nside - return npix2order(npix) + return cdshealpix.lonlat_to_healpix(ra, dec, order) -def ring2nest(nside: int, ipix: int) -> int: - order = nside2order(nside) +def ring2nest(order: int, ipix: int) -> int: return cdshealpix.from_ring(ipix, order) @@ -220,6 +210,5 @@ def order2mindist(order: np.ndarray | int) -> np.ndarray | float: np.ndarray of float or a single scalar float The minimum distance between pixels in arcminutes """ - pixel_nside = order2nside(order) - pixel_avgsize = nside2resol(pixel_nside, arcmin=True) + pixel_avgsize = order2resol(order, arcmin=True) return avgsize2mindist(pixel_avgsize) diff --git a/src/hats/pixel_math/partition_stats.py b/src/hats/pixel_math/partition_stats.py index 34f8c8f1..bc31345b 100644 --- a/src/hats/pixel_math/partition_stats.py +++ b/src/hats/pixel_math/partition_stats.py @@ -43,12 +43,10 @@ def generate_histogram( required_columns = [ra_column, dec_column] if not all(x in data.columns for x in required_columns): raise ValueError(f"Invalid column names in input: {ra_column}, {dec_column}") - mapped_pixels = hp.ang2pix( - 2**highest_order, + mapped_pixels = hp.radec2pix( + highest_order, data[ra_column].values, data[dec_column].values, - lonlat=True, - nest=True, ) mapped_pixel, count_at_pixel = np.unique(mapped_pixels, return_counts=True) histogram_result[mapped_pixel] += count_at_pixel.astype(np.int64) diff --git a/src/hats/pixel_math/spatial_index.py b/src/hats/pixel_math/spatial_index.py index 43c04b42..bc7ac6da 100644 --- a/src/hats/pixel_math/spatial_index.py +++ b/src/hats/pixel_math/spatial_index.py @@ -24,7 +24,7 @@ def compute_spatial_index(ra_values: List[float], dec_values: List[float]) -> np if len(ra_values) != len(dec_values): raise ValueError("ra and dec arrays should have the same length") - mapped_pixels = hp.ang2pix(2**SPATIAL_INDEX_ORDER, ra_values, dec_values, nest=True, lonlat=True) + mapped_pixels = hp.radec2pix(SPATIAL_INDEX_ORDER, ra_values, dec_values) return mapped_pixels diff --git a/tests/hats/pixel_math/test_healpix_shim.py b/tests/hats/pixel_math/test_healpix_shim.py index dd13654a..9d80e714 100644 --- a/tests/hats/pixel_math/test_healpix_shim.py +++ b/tests/hats/pixel_math/test_healpix_shim.py @@ -22,8 +22,7 @@ def test_mindist2avgsize2mindist(): def test_order2avgsize2order(): """Test that avgsize2order is inverse of hps.nside2resol(hps.order2nside, arcmin=True)""" order = np.arange(20) - nside = hps.order2nside(order) - assert_array_equal(hps.avgsize2order(hps.nside2resol(nside, arcmin=True)), order) + assert_array_equal(hps.avgsize2order(hps.order2resol(order, arcmin=True)), order) def test_margin2order(): @@ -97,92 +96,60 @@ def test_order2npix_invalid(): hps.order2npix(order) -def test_nside2pixarea(): +def test_order2pixarea(): orders = [0, 1, 5, 10, 20, 29] - nsides = [2**x for x in orders] npix = [12 * (4**order) for order in orders] pix_area_expected = [4 * np.pi / x for x in npix] - pix_area_test = [hps.nside2pixarea(nside) for nside in nsides] + pix_area_test = [hps.order2pixarea(order) for order in orders] assert pix_area_test == pix_area_expected -def test_nside2pixarea_degrees(): +def test_order2pixarea_degrees(): orders = [0, 1, 5, 10, 20, 29] - nsides = [2**x for x in orders] npix = [12 * (4**order) for order in orders] pix_area_expected = [np.rad2deg(np.rad2deg(4 * np.pi / x)) for x in npix] - pix_area_test = [hps.nside2pixarea(nside, degrees=True) for nside in nsides] + pix_area_test = [hps.order2pixarea(order, degrees=True) for order in orders] assert pix_area_test == pix_area_expected -def test_nside2resol(): +def test_order2resol(): orders = [0, 1, 5, 10, 20, 29] - nsides = [2**x for x in orders] - resol_expected = [np.sqrt(hps.nside2pixarea(nside)) for nside in nsides] - resol_test = [hps.nside2resol(nside) for nside in nsides] + resol_expected = [np.sqrt(hps.order2pixarea(order)) for order in orders] + resol_test = [hps.order2resol(order) for order in orders] assert resol_test == resol_expected -def test_nside2resol_arcmin(): +def test_order2resol_arcmin(): orders = [0, 1, 5, 10, 20, 29] nsides = [2**x for x in orders] - resol_expected = [np.rad2deg(np.sqrt(hps.nside2pixarea(nside))) * 60 for nside in nsides] - resol_test = [hps.nside2resol(nside, arcmin=True) for nside in nsides] + resol_expected = [np.rad2deg(np.sqrt(hps.order2pixarea(order))) * 60 for order in orders] + resol_test = [hps.order2resol(order, arcmin=True) for order in orders] assert resol_test == resol_expected -def test_ang2pix(): +def test_radec2pix_lonlat(): orders = [0, 1, 5, 10, 20, 29] ras = np.arange(-180.0, 180.0, 10.0) decs = np.arange(-90.0, 90.0, 180 // len(ras)) - colats = 90 - decs for order in orders: expected_pixels = cdshealpix.lonlat_to_healpix( Longitude(ras, unit="deg"), Latitude(decs, unit="deg"), order ) - pixels = hps.ang2pix(hps.order2nside(order), np.deg2rad(colats), np.deg2rad(ras), nest=True) + pixels = hps.radec2pix(order, ras, decs) assert np.all(pixels == expected_pixels) -def test_ang2pix_lonlat(): - orders = [0, 1, 5, 10, 20, 29] - ras = np.arange(-180.0, 180.0, 10.0) - decs = np.arange(-90.0, 90.0, 180 // len(ras)) - for order in orders: - expected_pixels = cdshealpix.lonlat_to_healpix( - Longitude(ras, unit="deg"), Latitude(decs, unit="deg"), order - ) - pixels = hps.ang2pix(hps.order2nside(order), ras, decs, nest=True, lonlat=True) - assert np.all(pixels == expected_pixels) - - -def test_ang2pix_ring(): - order = 4 - ras = np.arange(-180.0, 180.0, 10.0) - decs = np.arange(-90.0, 90.0, 180 // len(ras)) - colats = 90 - decs - with pytest.raises(NotImplementedError, match="RING"): - hps.ang2pix(hps.order2nside(order), ras, decs, lonlat=True) - with pytest.raises(NotImplementedError, match="RING"): - hps.ang2pix(hps.order2nside(order), colats, ras) - - -def test_ang2pix_invalid(): +def test_radec2pix_invalid(): orders = [0, 1, 5, 10, 20, 29] invalid_orders = [-1000, -1, 30, 40] ras = np.arange(-4000.0, 1000.0, 100.0) decs = np.arange(-1000.0, 1000.0, 2000.0 // len(ras)) - colats = 90 - decs for order in invalid_orders: with pytest.raises(ValueError, match="Invalid"): - hps.ang2pix(int(2**order), np.deg2rad(colats), np.deg2rad(ras), nest=True) - with pytest.raises(ValueError, match="Invalid"): - hps.ang2pix(int(2**order), ras, decs, nest=True, lonlat=True) + hps.radec2pix(order, ras, decs) for order in orders: with pytest.raises(ValueError, match="angle"): - hps.ang2pix(hps.order2nside(order), np.deg2rad(colats), np.deg2rad(ras), nest=True) - with pytest.raises(ValueError, match="angle"): - hps.ang2pix(hps.order2nside(order), ras, decs, nest=True, lonlat=True) + hps.radec2pix(order, ras, decs) def test_ring2nest(): @@ -190,5 +157,5 @@ def test_ring2nest(): for order in orders: ipix = np.arange(0, 12 * (4**order), 12 * (4**order) // 10) expected_ring_ipix = cdshealpix.from_ring(ipix, order) - test_ring_ipix = hps.ring2nest(2**order, ipix) + test_ring_ipix = hps.ring2nest(order, ipix) assert np.all(test_ring_ipix == expected_ring_ipix) diff --git a/tests/hats/pixel_math/test_spatial_index.py b/tests/hats/pixel_math/test_spatial_index.py index ca9b04eb..98377a07 100644 --- a/tests/hats/pixel_math/test_spatial_index.py +++ b/tests/hats/pixel_math/test_spatial_index.py @@ -16,12 +16,10 @@ def test_single(): """Single point. Adheres to specification.""" result = compute_spatial_index([5], [5]) - expected = hp.ang2pix( - 2**SPATIAL_INDEX_ORDER, + expected = hp.radec2pix( + SPATIAL_INDEX_ORDER, [5], [5], - nest=True, - lonlat=True, ) npt.assert_array_equal(result, expected) @@ -38,12 +36,10 @@ def test_short_list(): ra = [5, 1, 5] dec = [5, 1, 5] result = compute_spatial_index(ra, dec) - expected = hp.ang2pix( - 2**SPATIAL_INDEX_ORDER, + expected = hp.radec2pix( + SPATIAL_INDEX_ORDER, ra, dec, - nest=True, - lonlat=True, ) npt.assert_array_equal(result, expected) @@ -53,12 +49,10 @@ def test_list(): ra = [5, 5, 5, 1, 5, 5, 5, 1, 5] dec = [5, 5, 5, 1, 5, 5, 5, 1, 5] result = compute_spatial_index(ra, dec) - expected = hp.ang2pix( - 2**SPATIAL_INDEX_ORDER, + expected = hp.radec2pix( + SPATIAL_INDEX_ORDER, ra, dec, - nest=True, - lonlat=True, ) npt.assert_array_equal(result, expected) From f495c43a0577515dfec794f679ac87f1ac6aeaf0 Mon Sep 17 00:00:00 2001 From: Sean McGuire Date: Fri, 22 Nov 2024 11:36:30 -0500 Subject: [PATCH 15/19] remove unused variable --- tests/hats/pixel_math/test_healpix_shim.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/hats/pixel_math/test_healpix_shim.py b/tests/hats/pixel_math/test_healpix_shim.py index 9d80e714..7cc97f80 100644 --- a/tests/hats/pixel_math/test_healpix_shim.py +++ b/tests/hats/pixel_math/test_healpix_shim.py @@ -121,7 +121,6 @@ def test_order2resol(): def test_order2resol_arcmin(): orders = [0, 1, 5, 10, 20, 29] - nsides = [2**x for x in orders] resol_expected = [np.rad2deg(np.sqrt(hps.order2pixarea(order))) * 60 for order in orders] resol_test = [hps.order2resol(order, arcmin=True) for order in orders] assert resol_test == resol_expected From 48449119435a29defec1d7c1285ce8fd31e2ff44 Mon Sep 17 00:00:00 2001 From: Melissa DeLucchi <113376043+delucchi-cmu@users.noreply.github.com> Date: Fri, 22 Nov 2024 14:36:42 -0500 Subject: [PATCH 16/19] Return int64 from radec2pix (#432) --- src/hats/pixel_math/healpix_shim.py | 13 ++----------- 1 file changed, 2 insertions(+), 11 deletions(-) diff --git a/src/hats/pixel_math/healpix_shim.py b/src/hats/pixel_math/healpix_shim.py index 96eb18b7..2f1d8a22 100644 --- a/src/hats/pixel_math/healpix_shim.py +++ b/src/hats/pixel_math/healpix_shim.py @@ -58,14 +58,14 @@ def order2pixarea(order: int, degrees: bool = False) -> float: return pix_area_rad -def radec2pix(order: int, ra: float, dec: float) -> int: +def radec2pix(order: int, ra: float, dec: float) -> np.ndarray[np.int64]: if not is_order_valid(order): raise ValueError("Invalid value for order") ra = Longitude(ra, unit="deg") dec = Latitude(dec, unit="deg") - return cdshealpix.lonlat_to_healpix(ra, dec, order) + return cdshealpix.lonlat_to_healpix(ra, dec, order).astype(np.int64) def ring2nest(order: int, ipix: int) -> int: @@ -92,15 +92,6 @@ def ang2vec(ra, dec, **kwargs) -> np.ndarray: return np.array([coords.x.value, coords.y.value, coords.z.value]).T -## FITS -def read_map(*args, **kwargs): - return hp.read_map(*args, **kwargs) - - -def write_map(*args, **kwargs): - return hp.write_map(*args, **kwargs) - - ## Custom functions From e4ae032ed55787391569f340d27089535ae172f1 Mon Sep 17 00:00:00 2001 From: Sandro Campos Date: Fri, 22 Nov 2024 15:42:25 -0500 Subject: [PATCH 17/19] Migrate box filtering to use MOCpy (#428) * Use mocpy in box filtering * Update tests * Ensure moc is correct in box filter test * Remove unused healpy methods * Pin mocpy --- pyproject.toml | 2 +- .../healpix_dataset/healpix_dataset.py | 8 +- src/hats/pixel_math/box_filter.py | 98 ++--------- src/hats/pixel_math/healpix_shim.py | 16 -- src/hats/pixel_math/validators.py | 24 ++- tests/hats/catalog/test_catalog.py | 159 +++++------------- tests/hats/pixel_math/test_healpix_shim.py | 9 - 7 files changed, 68 insertions(+), 248 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/catalog/healpix_dataset/healpix_dataset.py b/src/hats/catalog/healpix_dataset/healpix_dataset.py index b5971717..f3d5f027 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 + 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 58c4b868..ad2d6a13 100644 --- a/src/hats/pixel_math/box_filter.py +++ b/src/hats/pixel_math/box_filter.py @@ -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 @@ -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: @@ -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)) diff --git a/src/hats/pixel_math/healpix_shim.py b/src/hats/pixel_math/healpix_shim.py index 2f1d8a22..ea2a903d 100644 --- a/src/hats/pixel_math/healpix_shim.py +++ b/src/hats/pixel_math/healpix_shim.py @@ -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 @@ -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 diff --git a/src/hats/pixel_math/validators.py b/src/hats/pixel_math/validators.py index 9f87d6f3..e81507ca 100644 --- a/src/hats/pixel_math/validators.py +++ b/src/hats/pixel_math/validators.py @@ -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) diff --git a/tests/hats/catalog/test_catalog.py b/tests/hats/catalog/test_catalog.py index 03883b2f..be96cdcc 100644 --- a/tests/hats/catalog/test_catalog.py +++ b/tests/hats/catalog/test_catalog.py @@ -6,6 +6,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 @@ -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 @@ -297,30 +297,32 @@ 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(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 + # Check that the previous filter is the same as intersecting the ra and dec filters 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()) + 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): # 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, 90)) filtered_pixels = filtered_catalog.get_healpix_pixels() assert len(filtered_pixels) == 2 @@ -332,13 +334,25 @@ 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 # 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, 90) def assert_is_subset_of(catalog, catalog_complement): pixels_catalog = catalog.get_healpix_pixels() @@ -346,138 +360,51 @@ 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 + # Some declination values are out of the [-90,90[ bounds 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 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)) - with pytest.raises(ValueError, match=ValidatorsErrors.INVALID_RADEC_RANGE): - small_sky_order1_catalog.filter_by_box(ra=(0, 360)) + # The declination values coincide 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 +521,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/pixel_math/test_healpix_shim.py b/tests/hats/pixel_math/test_healpix_shim.py index 7cc97f80..ba4b39ac 100644 --- a/tests/hats/pixel_math/test_healpix_shim.py +++ b/tests/hats/pixel_math/test_healpix_shim.py @@ -149,12 +149,3 @@ def test_radec2pix_invalid(): for order in orders: with pytest.raises(ValueError, match="angle"): hps.radec2pix(order, ras, decs) - - -def test_ring2nest(): - orders = [0, 1, 5, 10, 20, 29] - for order in orders: - ipix = np.arange(0, 12 * (4**order), 12 * (4**order) // 10) - expected_ring_ipix = cdshealpix.from_ring(ipix, order) - test_ring_ipix = hps.ring2nest(order, ipix) - assert np.all(test_ring_ipix == expected_ring_ipix) From fc48ad4a910092567b9e47f54b5ccfb0b1a61744 Mon Sep 17 00:00:00 2001 From: Melissa DeLucchi <113376043+delucchi-cmu@users.noreply.github.com> Date: Mon, 25 Nov 2024 14:48:42 -0500 Subject: [PATCH 18/19] Remove healpy references. (#433) --- docs/guide/contributing.rst | 14 +------------- pyproject.toml | 1 - src/hats/pixel_math/healpix_shim.py | 9 --------- src/hats/pixel_math/validators.py | 3 --- 4 files changed, 1 insertion(+), 26 deletions(-) diff --git a/docs/guide/contributing.rst b/docs/guide/contributing.rst index af3c0b34..b921fc06 100644 --- a/docs/guide/contributing.rst +++ b/docs/guide/contributing.rst @@ -58,19 +58,7 @@ Notes: .. tip:: Installing on Mac - - Native prebuilt binaries for healpy on Apple Silicon Macs - `do not yet exist `_, - so it's recommended to install via conda before proceeding to hats. - - .. code-block:: bash - - $ conda config --add channels conda-forge - $ conda install healpy - $ git clone https://github.com/astronomy-commons/hats - $ cd hats - $ pip install -e . - + When installing dev dependencies, make sure to include the single quotes. .. code-block:: bash diff --git a/pyproject.toml b/pyproject.toml index 7502b8b7..f8a66cf1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -22,7 +22,6 @@ dependencies = [ "aiohttp", # http filesystem support "astropy", "fsspec>=2023.10.0", # Used for abstract filesystems - "healpy", "jproperties", "matplotlib>=3.3,<3.9", "mocpy>=0.17.1", diff --git a/src/hats/pixel_math/healpix_shim.py b/src/hats/pixel_math/healpix_shim.py index ea2a903d..b5705a77 100644 --- a/src/hats/pixel_math/healpix_shim.py +++ b/src/hats/pixel_math/healpix_shim.py @@ -82,9 +82,6 @@ def ang2vec(ra, dec, **kwargs) -> np.ndarray: def avgsize2mindist(avg_size: np.ndarray) -> np.ndarray: """Get the minimum distance between pixels for a given average size - Average size is a "resolution" in healpy terms, i.e. given by - `healpy.nside2resol`. - We don't have the precise geometry of the healpix grid yet, so we are using average_size / mininimum_distance = 1.6 as a rough estimate. @@ -105,9 +102,6 @@ def avgsize2mindist(avg_size: np.ndarray) -> np.ndarray: def mindist2avgsize(mindist: np.ndarray) -> np.ndarray: """Get the average size for a given minimum distance between pixels - Average size is a "resolution" in healpy terms, i.e. given by - `healpy.nside2resol`. - We don't have the precise geometry of the healpix grid yet, so we are using average_size / mininimum_distance = 1.6 as a rough estimate. @@ -129,9 +123,6 @@ def mindist2avgsize(mindist: np.ndarray) -> np.ndarray: def avgsize2order(avg_size_arcmin: np.ndarray) -> np.ndarray: """Get the largest order with average healpix size larger than avg_size_arcmin - "Average" size is healpy's "resolution", so this function is - reverse to healpy.nside2resol(healpy.order2nside(order)).. - Parameters ---------- avg_size_arcmin : np.ndarray of float diff --git a/src/hats/pixel_math/validators.py b/src/hats/pixel_math/validators.py index e81507ca..f4928179 100644 --- a/src/hats/pixel_math/validators.py +++ b/src/hats/pixel_math/validators.py @@ -75,9 +75,6 @@ def validate_polygon(vertices: list[tuple[float, float]]): 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 From b14bc1467d02ca0187352fac8f6b88879d82091b Mon Sep 17 00:00:00 2001 From: Melissa DeLucchi Date: Mon, 25 Nov 2024 15:32:59 -0500 Subject: [PATCH 19/19] Address merge conflicts --- .github/workflows/pre-commit-ci.yml | 2 +- .github/workflows/testing-and-coverage.yml | 2 +- src/hats/pixel_math/validators.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/workflows/pre-commit-ci.yml b/.github/workflows/pre-commit-ci.yml index 642e7b9d..f31b3ee6 100644 --- a/.github/workflows/pre-commit-ci.yml +++ b/.github/workflows/pre-commit-ci.yml @@ -7,7 +7,7 @@ on: push: branches: [ main ] pull_request: - branches: [ main, margin ] + branches: [ main ] jobs: pre-commit-ci: diff --git a/.github/workflows/testing-and-coverage.yml b/.github/workflows/testing-and-coverage.yml index 8438b944..3d3b867f 100644 --- a/.github/workflows/testing-and-coverage.yml +++ b/.github/workflows/testing-and-coverage.yml @@ -7,7 +7,7 @@ on: push: branches: [ main ] pull_request: - branches: [ main, margin ] + branches: [ main ] jobs: build: diff --git a/src/hats/pixel_math/validators.py b/src/hats/pixel_math/validators.py index def08b2e..9f61601d 100644 --- a/src/hats/pixel_math/validators.py +++ b/src/hats/pixel_math/validators.py @@ -78,7 +78,7 @@ def check_polygon_is_valid(vertices: np.ndarray): Args: vertices (np.ndarray): The polygon vertices, in cartesian coordinates """ - vertices_xyz = hp.ang2vec(*vertices.T, lonlat=True) + vertices_xyz = hp.ang2vec(*vertices.T) # Compute the normal between each pair of neighboring vertices second_vertices = np.roll(vertices_xyz, -1, axis=0)