Implement range to moc function

healpix_alchemy/tests/benchmarks/conftest.py

@@ -40,3 +40,8 @@ def random_fields(cursor, tables, request):
 @pytest.fixture
 def random_sky_map(cursor, tables):
     return data.get_random_sky_map(20_000, cursor)
+
+
+@pytest.fixture
+def random_moc_from_cone(cursor, tables):
+    return data.get_random_moc_from_cone(15, cursor)

healpix_alchemy/tests/benchmarks/data.py

@@ -7,21 +7,23 @@
 """
 import io
 
-from astropy.coordinates import SkyCoord, uniform_spherical_random_surface
+from astropy.coordinates import (SkyCoord, uniform_spherical_random_surface,
+                                 Angle, Longitude, Latitude)
 from astropy import units as u
 from mocpy import MOC
 import numpy as np
 import pytest
 
 from ...constants import HPX, LEVEL, PIXEL_AREA
 from ...types import Tile
-from .models import Galaxy, Field, FieldTile, Skymap, SkymapTile
+from .models import Galaxy, Field, FieldTile, Skymap, SkymapTile, TileList
 
 (
     RANDOM_GALAXIES_SEED,
     RANDOM_FIELDS_SEED,
-    RANDOM_SKY_MAP_SEED
-) = np.random.SeedSequence(12345).spawn(3)
+    RANDOM_SKY_MAP_SEED,
+    RANDOM_MOC_SEED
+) = np.random.SeedSequence(12345).spawn(4)
 
 
 def get_ztf_footprint_corners():
@@ -122,3 +124,17 @@ def get_random_sky_map(n, cursor):
     cursor.copy_from(f, SkymapTile.__tablename__)
 
     return tiles, probdensity
+
+
+def get_random_moc_from_cone(depth, cursor):
+    rng = np.random.default_rng(RANDOM_MOC_SEED)
+    # Randomly choose lon, lat & radius for cone search
+    lon = Longitude(360 * rng.random() * u.deg)
+    lat = Latitude((180 * rng.random() - 90) * u.deg)
+    radius = Angle(rng.normal(loc=10, scale=2.5) * u.deg)
+    moc = MOC.from_cone(lon, lat, radius, depth)
+
+    f = io.StringIO('\n'.join(Tile.tiles_from_moc(moc)))
+    cursor.copy_from(f, TileList.__tablename__)
+
+    return moc

healpix_alchemy/tests/benchmarks/models.py

@@ -35,3 +35,7 @@ class SkymapTile(Base):
     id = sa.Column(sa.ForeignKey(Skymap.id), primary_key=True)
     hpx = sa.Column(Tile, primary_key=True, index=True)
     probdensity = sa.Column(sa.Float, nullable=False)
+
+
+class TileList(Base):
+    hpx = sa.Column(Tile, primary_key=True, index=True)

healpix_alchemy/tests/benchmarks/test_benchmarks.py

@@ -3,9 +3,11 @@
 import numpy as np
 import sqlalchemy as sa
 import pytest
+import re
 
 from ... import func
-from .models import Galaxy, FieldTile, SkymapTile
+from .models import Galaxy, FieldTile, SkymapTile, TileList
+from ...types import Tile
 
 
 @pytest.fixture
@@ -104,3 +106,17 @@ def test_fields_in_90pct_credible_region(bench, random_fields, random_sky_map):
 
     # Run benchmark
     bench(query)
+
+
+def test_to_moc(bench, random_moc_from_cone):
+    # Assemble Query
+    query = sa.select(TileList.hpx)
+    rangeset = map(lambda x:
+                   (int(re.split("[,]|[(]", str(x))[2]),
+                    int(re.split("[,]|[(]", str(x))[3])),
+                   bench(query))
+
+    # Expected result
+    expected = random_moc_from_cone
+
+    assert expected == Tile.moc_from_tiles(rangeset, 2**29)

healpix_alchemy/types.py

@@ -80,6 +80,13 @@ def tiles_from_polygon_skycoord(cls, polygon):
     def tiles_from_moc(cls, moc):
         return (f'[{lo},{hi})' for lo, hi in moc._interval_set.nested)
 
+    @classmethod
+    def moc_from_tiles(cls, rangeSet, nside):
+        depth = int(np.log2(nside))
+        MOCstr = f'{depth}/' + ' '.join(map(lambda x: f'{x[0]}-{x[1]-1}',
+                                            rangeSet))
+        return MOC.from_string(MOCstr)
+
 
 @sa.event.listens_for(sa.Index, 'after_parent_attach')
 def _create_indices(index, parent):