Implement range to moc function

healpix_alchemy/tests/test_types.py

@@ -0,0 +1,57 @@
+import numpy as np
+from mocpy import MOC
+import itertools
+from astropy_healpix import HEALPix
+
+from .. import types
+
+
+def random_MOC(max_l):
+    json = {str(level): list(
+        np.unique(
+            np.random.randint(0, 12 * (4 ** level),
+                              np.random.randint(12 * (4 ** level)))))
+            for level in range(max_l)}
+    return MOC.from_json(json), json
+
+
+def to_ranges(iterable):
+    iterable = sorted(set(iterable))
+    for key, group in itertools.groupby(enumerate(iterable),
+                                        lambda t: t[1] - t[0]):
+        group = list(group)
+        yield group[0][1], group[-1][1]+1
+
+
+def test_to_moc_ring():
+    l_lim = 12
+    moc, json = random_MOC(l_lim)
+    for level in range(l_lim-1):
+        for pixel in json[str(level)]:
+            for i in range(4):
+                json[str(level + 1)].append(pixel * 4 + i)
+        json[str(level + 1)] = np.unique(json[str(level + 1)])
+    nested_hpx_list = json[str(l_lim-1)]
+    hp = HEALPix(nside=2**(l_lim-1), order='ring')
+    ring_hpx_list = hp.nested_to_ring(nested_hpx_list)
+    ring_hpx_ranges = to_ranges(ring_hpx_list)
+    point = types.Point()
+    assert point.to_moc(rangeSet=ring_hpx_ranges,
+                        nside=2**(l_lim-1),
+                        index='ring') == moc
+
+
+def test_to_moc_nested():
+    l_lim = 12
+    moc, json = random_MOC(l_lim)
+    for level in range(l_lim-1):
+        for pixel in json[str(level)]:
+            for i in range(4):
+                json[str(level + 1)].append(pixel * 4 + i)
+        json[str(level + 1)] = np.unique(json[str(level + 1)])
+    nested_hpx_list = json[str(l_lim-1)]
+    nested_hpx_ranges = to_ranges(nested_hpx_list)
+    point = types.Point()
+    assert point.to_moc(rangeSet=nested_hpx_ranges,
+                        nside=2**(l_lim-1),
+                        index='nested') == moc

healpix_alchemy/types.py

@@ -3,7 +3,7 @@
 from numbers import Integral
 
 from astropy.coordinates import SkyCoord
-from astropy_healpix import uniq_to_level_ipix
+from astropy_healpix import uniq_to_level_ipix, HEALPix
 from mocpy import MOC
 import numpy as np
 import sqlalchemy as sa
@@ -28,6 +28,15 @@ def process_bind_param(self, value, dialect):
             value = int(value)
         return value
 
+    def to_moc(self, rangeSet, nside, index):
+        healpixSet = np.unique(np.concatenate([np.arange(start, stop, 1) for
+                                               (start, stop) in rangeSet]))
+        if index.lower() == "ring":
+            hp = HEALPix(nside=nside, order='ring')
+            healpixSet = hp.ring_to_nested(healpixSet)
+        depth = np.ones(np.shape(healpixSet)) * np.log2(nside)
+        return MOC.from_healpix_cells(healpixSet, depth)
+
 
 class Tile(sa.TypeDecorator):