Add support for pixel ranges and composable geometries.

.github/workflows/python-package-pr.yml

@@ -24,9 +24,7 @@ jobs:
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
-        pip install flake8 pytest
-        if [ -f requirements.txt ]; then pip install -r requirements.txt; fi
-        pip install .[full,parquet,healpy]
+        pip install .[parquet,healpy,test]
     - name: Lint with flake8
       run: |
         flake8

.github/workflows/python-package-publish.yml

@@ -25,9 +25,7 @@ jobs:
       - name: Build and install
         run: |
           python -m pip install --upgrade pip setuptools
-          python -m pip install pytest
-          if [ -f requirements.txt ]; then pip install -r requirements.txt; fi
-          python -m pip install .[full,parquet,healpy]
+          python -m pip install .[parquet,healpy,test]
       - name: Run tests
         run: |
           cd tests

.github/workflows/python-package-push.yml

@@ -24,9 +24,7 @@ jobs:
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
-        pip install flake8 pytest
-        if [ -f requirements.txt ]; then pip install -r requirements.txt; fi
-        pip install .[full,parquet,healpy]
+        python -m pip install .[parquet,healpy,test]
     - name: Lint with flake8
       run: |
         flake8

.readthedocs.yml

@@ -5,19 +5,19 @@
 # Required
 version: 2
 
+# Set the version of Python and other tools you might need
+build:
+  os: ubuntu-22.04
+  tools:
+    python: "3.11"
+
 # Build documentation in the docs/ directory with Sphinx
 sphinx:
   builder: html
   configuration: docs/conf.py
   fail_on_warning: false
 
-# Optionally build your docs in additional formats such as PDF
-formats:
-  - pdf
-
-# Optionally set the version of Python and requirements required to build your docs
 python:
-  version: 3.7
   install:
     - requirements: docs/requirements.txt
     - requirements: requirements.txt

docs/basic_interface.rst

@@ -308,3 +308,5 @@ For example, we can take render our map as a collection of hexagonal cells using
     plt.hexbin(ra, dec, C=hsp_map[vpix])
     plt.colorbar()
     plt.show()
+
+For more sophisticated visualizations and projections of HealSparse maps, please see `SkyProj <https://skyproj.readthedocs.io/en/latest/>`_.

docs/geometry.rst

@@ -4,9 +4,12 @@
 HealSparse Geometry
 ===================
 
-:code:`HealSparse` has a basic geometry library that allows you to generate maps from circles and convex polygons, as supported by :code:`healpy`.  Each geometric object is associated with a single value.  On construction, geometry objects only contain information about the shape, and they are only rendered onto a `HEALPix` grid when requested.
+:code:`HealSparse` has a basic geometry library that allows you to generate maps from circles, convex polygons, and ellipses as supported by :code:`hpgeom`.  Each geometric object is associated with a single value.  On construction, geometry objects only contain information about the shape, and they are only rendered onto a `HEALPix` grid when requested.
 
-There are two methods to realize geometry objects.  The first is that each object can be used to generate a :code:`HealSparseMap` map, and the second, for integer-valued objects is the :code:`realize_geom()` method which can be used to combine multiple objects by :code:`or`-ing the integer values together.
+There are a few methods to realize geometry objects.
+The easiest is to combine a geometric object with a :code:`HealSparseMap` map, with the ``or``, ``and``, or ``add`` operation.
+One can generate a :code:`HealSparseMap` from the geometric object.
+Finally, for integer-value objects one can use the :code:`realize_geom()` method to combine multiple objects by :cod:`or`-ing the integer values together.
 
 
 HealSparse Geometry Shapes
@@ -46,6 +49,30 @@ The three shapes supported are :code:`Circle`, :code:`Ellipse`, and :code:`Polyg
                               value=8)
 
 
+Combining Geometric Objects with Maps
+-------------------------------------
+
+Given a map, it is very simple to combine geometric objects to build up complex shapes/masks/etc.
+Behind the scenes, large geometric objects are rendered with :code:`hpgeom` pixel ranges which leads to greater memory efficiency.
+
+.. code-block :: python
+
+    import healsparse
+    import numpy as np
+
+    # Create an empty map.
+    m = healsparse.HealSparseMap.make_empty(32, 4096, np.uint16)
+
+    # Set a large circle to a value using the ``or`` operation
+    m |= healsparse.Circle(ra=200.0, dec=20.0, radius=5.0, value=1)
+
+    # Remove a small circle from the center using the ``and`` operation
+    m &= healsparse.Circle(ra=200.0, dec=20.0, radius=1.0, value=0)
+
+    # And add in another circle.
+    m += healsparse.Circle(ra=202.0, dec=21.0, radius=0.5, value=10)
+
+
 Making a Map
 ------------
 

docs/install.rst

@@ -9,7 +9,7 @@ Installation of `healpy <https://github.com/healpy/healpy>`_ is required for rea
 
 .. code-block:: bash
 
-  conda install -c conda-forge hpgeom
+  conda install -c conda-forge healsparse
 
 or
 

healsparse/fits_shim.py

@@ -241,7 +241,7 @@ def _write_filename(filename, c_hdr, s_hdr, cov_index_map, sparse_map,
 
     integer_map = sparse_map.dtype.fields is None and is_integer_value(sparse_map.dtype.type(0))
     if integer_map:
-        compression = "RICE"
+        compression = "RICE_1"
     else:
         compression = "GZIP_2"
 

healsparse/geom.py

@@ -1,6 +1,5 @@
 import numpy as np
 import hpgeom as hpg
-from .healSparseMap import HealSparseMap
 from .utils import is_integer_value
 import numbers
 
@@ -16,21 +15,19 @@ def realize_geom(geom, smap, type='or'):
     smap : `HealSparseMap`
         The map in which to realize the objects.
     type : `str`
-        Way to combine the list of geometric objects.  Default
-        is to "or" them.
+        Way to combine the list of geometric objects.
+        Currently only supports ``or``.
     """
-
-    if type != 'or':
-        raise ValueError('type of composition must be or')
+    if type not in ['or']:
+        raise ValueError('Type of composition must be ``or``')
 
     if not smap.is_integer_map:
-        raise ValueError('can only or geometry objects into an integer map')
+        raise ValueError(f'Can only {type} geometry objects into an integer map')
 
     if not isinstance(geom, (list, tuple)):
         geom = [geom]
 
-    # split the geom objects up by value
-    gdict = {}
+    # Check all the values before starting.
     for g in geom:
         value = g.value
         if isinstance(value, (tuple, list, np.ndarray)):
@@ -44,35 +41,9 @@ def realize_geom(geom, smap, type='or'):
             _check_int(value)
             _check_int_size(value, smap.dtype)
 
-        if value not in gdict:
-            gdict[value] = [g]
-        else:
-            gdict[value].append(g)
-
-    # deal with each value separately and add to
-    # the map
-    for value, glist in gdict.items():
-        for i, g in enumerate(glist):
-            tpixels = g.get_pixels(nside=smap.nside_sparse)
-            if i == 0:
-                pixels = tpixels
-            else:
-                oldsize = pixels.size
-                newsize = oldsize + tpixels.size
-                # need refcheck=False because it will fail when running
-                # the python profiler; I infer that the profiler holds
-                # a reference to the objects
-                pixels.resize(newsize, refcheck=False)
-                pixels[oldsize:] = tpixels
-
-        pixels = np.unique(pixels)
-
-        if smap.is_wide_mask_map:
-            smap.set_bits_pix(pixels, value)
-        else:
-            values = smap.get_values_pix(pixels)
-            values |= value
-            smap.update_values_pix(pixels, values)
+    # Now generate the map.
+    for g in geom:
+        smap |= g
 
 
 def _check_int(x):
@@ -120,6 +91,17 @@ def get_pixels(self, *, nside):
         """
         raise NotImplementedError('Implement get_pixels')
 
+    def get_pixel_ranges(self, *, nside):
+        """
+        Get pixel ranges for this geometric shape.
+
+        Parameters
+        ----------
+        nside : `int`
+            HEALPix nside for the pixels.
+        """
+        raise NotImplementedError("Implement get_pixel_ranges")
+
     def get_map(self, *, nside_coverage, nside_sparse, dtype, wide_mask_maxbits=None):
         """
         Get a healsparse map corresponding to this geometric primitive.
@@ -139,6 +121,7 @@ def get_map(self, *, nside_coverage, nside_sparse, dtype, wide_mask_maxbits=None
         -------
         hsmap : `healsparse.HealSparseMap`
         """
+        from .healSparseMap import HealSparseMap
 
         x = np.zeros(1, dtype=dtype)
         if is_integer_value(x[0]):
@@ -187,6 +170,7 @@ def get_map_like(self, sparseMap):
         -------
         hsmap : `healsparse.HealSparseMap`
         """
+        from .healSparseMap import HealSparseMap
 
         if not isinstance(sparseMap, HealSparseMap):
             raise RuntimeError("Input sparseMap must be a HealSparseMap")
@@ -258,6 +242,17 @@ def get_pixels(self, *, nside):
             inclusive=False,
         )
 
+    def get_pixel_ranges(self, *, nside):
+        return hpg.query_circle(
+            nside,
+            self._ra,
+            self._dec,
+            self._radius,
+            nest=True,
+            inclusive=False,
+            return_pixel_ranges=True,
+        )
+
     def __repr__(self):
         s = 'Circle(ra=%.16g, dec=%.16g, radius=%.16g, value=%s)'
         return s % (self._ra, self._dec, self._radius, repr(self._value))
@@ -315,15 +310,23 @@ def vertices(self):
         return self._vertices
 
     def get_pixels(self, *, nside):
-        pixels = hpg.query_polygon(
+        return hpg.query_polygon(
             nside,
             self._ra,
             self._dec,
             nest=True,
             inclusive=False,
         )
 
-        return pixels
+    def get_pixel_ranges(self, *, nside):
+        return hpg.query_polygon(
+            nside,
+            self._ra,
+            self._dec,
+            nest=True,
+            inclusive=False,
+            return_pixel_ranges=True,
+        )
 
     def __repr__(self):
         ras = repr(self._ra)
@@ -413,7 +416,20 @@ def get_pixels(self, *, nside):
             self._semi_minor,
             self._alpha,
             nest=True,
-            inclusive=False
+            inclusive=False,
+        )
+
+    def get_pixel_ranges(self, *, nside):
+        return hpg.query_ellipse(
+            nside,
+            self._ra,
+            self._dec,
+            self._semi_major,
+            self._semi_minor,
+            self._alpha,
+            nest=True,
+            inclusive=False,
+            return_pixel_ranges=True,
         )
 
     def __repr__(self):