diff --git a/src/hipscat/pixel_math/healpix_shim.py b/src/hipscat/pixel_math/healpix_shim.py index 381956d7..4b4a34e4 100644 --- a/src/hipscat/pixel_math/healpix_shim.py +++ b/src/hipscat/pixel_math/healpix_shim.py @@ -1,4 +1,5 @@ import healpy as hp +import numpy as np # pylint: disable=missing-function-docstring @@ -98,3 +99,95 @@ def read_map(*args, **kwargs): def write_map(*args, **kwargs): return hp.write_map(*args, **kwargs) + + +## Custom functions + + +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. + + Parameters + ---------- + avg_size : np.ndarray of float + The average size of a healpix pixel + + Returns + ------- + np.ndarray of float + The minimum distance between pixels for the given average size + """ + return avg_size / 1.6 + + +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. + + Parameters + ---------- + mindist : np.ndarray of float + The minimum distance between pixels + + Returns + ------- + np.ndarray of float + The average size of a healpix pixel for the given minimum distance + between pixels. + """ + return mindist * 1.6 + + +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 + The average size of a healpix pixel in arcminutes + + Returns + ------- + np.ndarray of int + The largest healpix order for which the average size is larger than avg_size_arcmin + """ + # resolution = sqrt(4pi / (12 * 4**order)) => order = log2(sqrt(4pi / 12) / resolution) + order_float = np.log2(np.sqrt(np.pi / 3) / np.radians(avg_size_arcmin / 60.0)) + return np.array(order_float, dtype=int) + + +def margin2order(margin_thr_arcmin: np.ndarray) -> np.ndarray: + """Get the largest order for which distance between pixels is less than margin_thr_arcmin + + We don't have the precise geometry of the healpix grid yet, + we are using average_size / mininimum_distance = 1.6 + as a rough estimate. + + Parameters + ---------- + margin_thr_arcmin : np.ndarray of float + The minimum distance between pixels in arcminutes + + Returns + ------- + np.ndarray of int + The largest healpix order for which the distance between pixels is less than margin_thr_arcmin + """ + avg_size_arcmin = mindist2avgsize(margin_thr_arcmin) + return avgsize2order(avg_size_arcmin) diff --git a/tests/hipscat/pixel_math/test_healpix_shim.py b/tests/hipscat/pixel_math/test_healpix_shim.py new file mode 100644 index 00000000..68f626e8 --- /dev/null +++ b/tests/hipscat/pixel_math/test_healpix_shim.py @@ -0,0 +1,30 @@ +import numpy as np +from numpy.testing import assert_allclose, assert_array_equal + +from hipscat.pixel_math import healpix_shim as hps + + +def test_avgsize2mindist2avgsize(): + """Test that avgsize2mindist is inverse of mindist2avgsize""" + avgsize = np.logspace(-5, 5, 21) + assert_allclose(hps.mindist2avgsize(hps.avgsize2mindist(avgsize)), avgsize, rtol=1e-12) + + +def test_mindist2avgsize2mindist(): + """Test that mindist2avgsize is inverse of avgsize2mindist""" + mindist = np.logspace(-3.2, 2.8, 21) + assert_allclose(hps.avgsize2mindist(hps.mindist2avgsize(mindist)), mindist, rtol=1e-12) + + +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) + + +def test_margin2order(): + """Test margin2order for some pre-computed values""" + margin_thr_arcmin = np.array([1 / 60, 10 / 60, 1, 5, 60]) + orders = np.array([17, 13, 11, 8, 5]) + assert_array_equal(hps.margin2order(margin_thr_arcmin), orders)