Merge pull request #318 from astronomy-commons/avgsize_to_order

Add utility healpix functions to convert between order, average pixel size and minimum distance between edges

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)

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)