Merge pull request #159 from LSSTDESC/geom_nside_render

Add nside_render to geom objects to allow for rendering at a different (coarser) resolution than the corresponding map.

docs/geometry.rst

@@ -4,7 +4,15 @@
 HealSparse Geometry
 ===================
 
-: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.
+:code:`HealSparse` has a basic geometry library that allows you to generate maps from circles, convex polygons, ellipses, and boxes 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.
+
+In addition to information about the shape itself, a geometric object may optionally contain a value of :code:`nside_render`.
+This indicates that the shape should always be rendered at this given resolution, no matter the resolution of the map that it is being combined with.
+(Note that you can only render at a resolution that is less than or equal to the map resolution, or else a :code:`ValueError` is raised.)
+This functionality may be useful if one is building a map that may be used with multiple resolutions, and one wants to ensure that a higher and lower resolution maps have exactly the same outline for these shapes.
+If no :code:`nside_render` is set with the object it will always be rendered at the same resolution as the corresponding map or via the :code:`nside` parameter of :code:`get_pixels()` and :code:`get_pixel_ranges()`.
 
 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.
@@ -15,7 +23,8 @@ Finally, for integer-value objects one can use the :code:`realize_geom()` method
 HealSparse Geometry Shapes
 --------------------------
 
-The three shapes supported are :code:`Circle`, :code:`Ellipse`, and :code:`Polygon`.  They share a base class, and while the instantiation is different, the operations are the same.
+The four shapes supported are :code:`Circle`, :code:`Ellipse`, :code:`Polygon`, and :code:`Box`.
+They share a base class, and while the instantiation is different, the operations are the same.
 
 **Circle**
 
@@ -49,11 +58,23 @@ The three shapes supported are :code:`Circle`, :code:`Ellipse`, and :code:`Polyg
                               value=8)
 
 
+**Box**
+
+.. code-block :: python
+
+    # All units are decimal degrees
+    # A box differs from a polygon in that the sides will be constant ra/dec rather
+    # than great circles.
+    # See https://hpgeom.readthedocs.io/en/latest .
+    box = healsparse.Box(ra1=20.0, ra2=30.0, dec1=10.0, dec2=5.0, value=True)
+
+
 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.
+Note that these operations can be applied to integer or boolean maps.
 
 .. code-block :: python
 
@@ -89,7 +110,8 @@ To make a map from a geometry object, use the :code:`get_map()` method as such.
 Using :code:`realize_geom()`
 ----------------------------
 
-You can only use :code:`realize_geom()` to create maps from combinations of polygons if you are using integer maps, and want to :code:`or` them together.  This method is more memory efficient than generating each individual individual map and combining them, as above.
+You can only use :code:`realize_geom()` to create maps from combinations of polygons if you are using integer maps, and want to :code:`or` them together.
+This method is more memory efficient than generating each individual individual map and combining them, as above.
 
 .. code-block :: python
 

healsparse/__init__.py

@@ -21,6 +21,6 @@
 from .operations import ufunc_union, ufunc_intersection
 from .operations import divide_intersection, floor_divide_intersection
 from . import geom
-from .geom import Circle, Ellipse, Polygon, realize_geom
+from .geom import Box, Circle, Ellipse, Polygon, realize_geom
 from .utils import WIDE_MASK
 from .cat_healsparse_files import cat_healsparse_files

healsparse/geom.py

@@ -88,8 +88,25 @@ def get_pixels(self, *, nside):
         ----------
         nside : `int`
             HEALPix nside for the pixels.
+
+        Raises
+        ------
+        ValueError : If shape has nside_render set, this is raised if
+            nside < nside_render.
         """
-        raise NotImplementedError('Implement get_pixels')
+        if self._nside_render is not None:
+            if nside < self._nside_render:
+                raise ValueError(f"Cannot render a Circle with {self._nside_render} into nside={nside}")
+            _nside = self._nside_render
+        else:
+            _nside = nside
+
+        pixels = self._render(nside_render=_nside, return_pixel_ranges=False)
+
+        if self._nside_render is not None:
+            return hpg.upgrade_pixels(_nside, pixels, nside)
+        else:
+            return pixels
 
     def get_pixel_ranges(self, *, nside):
         """
@@ -99,8 +116,41 @@ def get_pixel_ranges(self, *, nside):
         ----------
         nside : `int`
             HEALPix nside for the pixels.
+
+        Raises
+        ------
+        ValueError : If shape has nside_render set, this is raised if
+            nside < nside_render.
+        """
+        if self._nside_render is not None:
+            if nside < self._nside_render:
+                raise ValueError(f"Cannot render a Circle with {self._nside_render} into nside={nside}")
+            _nside = self._nside_render
+        else:
+            _nside = nside
+
+        pixel_ranges = self._render(nside_render=_nside, return_pixel_ranges=True)
+
+        if self._nside_render is not None:
+            return hpg.upgrade_pixel_ranges(_nside, pixel_ranges, nside)
+        else:
+            return pixel_ranges
+
+    def _render(self, *, nside_render, return_pixel_ranges):
+        """Internal method to render to pixels/ranges for this shape.
+
+        Parameters
+        ----------
+        nside_render : `int`
+            Rendering resolution.
+        return_pixel_ranges : `bool`
+            Return pixel ranges instead of pixels?
+
+        Returns
+        -------
+        pixels or pixel_ranges : `np.ndarray`
         """
-        raise NotImplementedError("Implement get_pixel_ranges")
+        raise NotImplementedError("The _render method must be overridden.")
 
     def get_map(self, *, nside_coverage, nside_sparse, dtype, wide_mask_maxbits=None):
         """
@@ -126,6 +176,8 @@ def get_map(self, *, nside_coverage, nside_sparse, dtype, wide_mask_maxbits=None
         x = np.zeros(1, dtype=dtype)
         if is_integer_value(x[0]):
             sentinel = 0
+        elif dtype == np.bool_ or dtype == bool:
+            sentinel = False
         else:
             sentinel = hpg.UNSEEN
 
@@ -199,12 +251,17 @@ class Circle(GeomBase):
         Radius in degrees (scalar-only).
     value : number
         Value for pixels in the map (scalar or list of bits for `wide_mask`)
+    nside_render : `int`, optional
+        If this is set, the shape will always be rendered at this
+        nside and then these pixels will be 'upgraded' to the resolution
+        of the map.
     """
-    def __init__(self, *, ra, dec, radius, value):
+    def __init__(self, *, ra, dec, radius, value, nside_render=None):
         self._ra = ra
         self._dec = dec
         self._radius = radius
         self._value = value
+        self._nside_render = nside_render
         sc_ra = np.isscalar(self._ra)
         sc_dec = np.isscalar(self._dec)
         sc_radius = np.isscalar(self._radius)
@@ -232,30 +289,20 @@ def radius(self):
         """
         return self._radius
 
-    def get_pixels(self, *, nside):
+    def _render(self, *, nside_render, return_pixel_ranges):
         return hpg.query_circle(
-            nside,
+            nside_render,
             self._ra,
             self._dec,
             self._radius,
             nest=True,
             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,
+            return_pixel_ranges=return_pixel_ranges,
         )
 
     def __repr__(self):
-        s = 'Circle(ra=%.16g, dec=%.16g, radius=%.16g, value=%s)'
-        return s % (self._ra, self._dec, self._radius, repr(self._value))
+        s = 'Circle(ra=%.16g, dec=%.16g, radius=%.16g, value=%s, nside_render=%s)'
+        return s % (self._ra, self._dec, self._radius, repr(self._value), repr(self._nside_render))
 
 
 class Polygon(GeomBase):
@@ -273,7 +320,7 @@ class Polygon(GeomBase):
     value : number
         Value for pixels in the map
     """
-    def __init__(self, *, ra, dec, value):
+    def __init__(self, *, ra, dec, value, nside_render=None):
         ra = np.array(ra, ndmin=1)
         dec = np.array(dec, ndmin=1)
 
@@ -285,6 +332,7 @@ def __init__(self, *, ra, dec, value):
         self._dec = dec
         self._vertices = hpg.angle_to_vector(ra, dec, lonlat=True)
         self._value = value
+        self._nside_render = nside_render
 
         self._is_integer = is_integer_value(value)
 
@@ -309,31 +357,22 @@ def vertices(self):
         """
         return self._vertices
 
-    def get_pixels(self, *, nside):
+    def _render(self, *, nside_render, return_pixel_ranges):
         return hpg.query_polygon(
-            nside,
+            nside_render,
             self._ra,
             self._dec,
             nest=True,
             inclusive=False,
-        )
-
-    def get_pixel_ranges(self, *, nside):
-        return hpg.query_polygon(
-            nside,
-            self._ra,
-            self._dec,
-            nest=True,
-            inclusive=False,
-            return_pixel_ranges=True,
+            return_pixel_ranges=return_pixel_ranges,
         )
 
     def __repr__(self):
         ras = repr(self._ra)
         decs = repr(self._dec)
 
-        s = 'Polygon(ra=%s, dec=%s, value=%s)'
-        return s % (ras, decs, repr(self._value))
+        s = 'Polygon(ra=%s, dec=%s, value=%s, nside_render=%s)'
+        return s % (ras, decs, repr(self._value), repr(self._nside_render))
 
 
 class Ellipse(GeomBase):
@@ -355,13 +394,14 @@ class Ellipse(GeomBase):
     value : number
         Value for pixels in the map (scalar or list of bits for `wide_mask`).
     """
-    def __init__(self, *, ra, dec, semi_major, semi_minor, alpha, value):
+    def __init__(self, *, ra, dec, semi_major, semi_minor, alpha, value, nside_render=None):
         self._ra = ra
         self._dec = dec
         self._semi_major = semi_major
         self._semi_minor = semi_minor
         self._alpha = alpha
         self._value = value
+        self._nside_render = nside_render
         sc_ra = np.isscalar(self._ra)
         sc_dec = np.isscalar(self._dec)
         sc_semi_major = np.isscalar(self._semi_major)
@@ -407,31 +447,95 @@ def alpha(self):
         """
         return self._alpha
 
-    def get_pixels(self, *, nside):
+    def _render(self, *, nside_render, return_pixel_ranges):
         return hpg.query_ellipse(
-            nside,
+            nside_render,
             self._ra,
             self._dec,
             self._semi_major,
             self._semi_minor,
             self._alpha,
             nest=True,
             inclusive=False,
+            return_pixel_ranges=return_pixel_ranges,
         )
 
-    def get_pixel_ranges(self, *, nside):
-        return hpg.query_ellipse(
-            nside,
-            self._ra,
-            self._dec,
-            self._semi_major,
-            self._semi_minor,
-            self._alpha,
+    def __repr__(self):
+        s = ("Ellipse(ra=%.16g, dec=%16g, semi_major=%16g, semi_minor=%16g, alpha=%16g, value=%s, "
+             "nside_render=%s)")
+        return s % (self._ra, self._dec, self._semi_major, self._semi_minor, self._alpha, repr(self._value),
+                    repr(self._nside_render))
+
+
+class Box(GeomBase):
+    """
+    A geometric shape that has sides of constant lon/lat.
+
+    This shape is in contrast with a Polygon which will have great
+    circle boundaries. See hpgeom.query_box() for details.
+
+    Parameters
+    ----------
+    ra1, ra2 : `float`
+        RA in degrees. All points within [ra1, ra2] will be selected.
+        If ra1 > ra2 then the box will wrap around 360 degrees. If
+        ra1 == 0.0 and ra2 == 360.0 then the box will contain points at
+        all right ascensions.
+    dec1, dec2 : `float`
+        Declination in degrees. All points within [dec1, dec2] will be
+        selected.  If dec1 or dec2 is 90.0 or -90.0 then the box will
+        be an arc of a circle with the center at the north/south pole.
+    value : number
+        Value for pixels in the map (scalar or list of bits for `wide_mask`)
+    nside_render : `int`, optional
+        If this is set, the shape will always be rendered at this
+        nside and then these pixels will be 'upgraded' to the resolution
+        of the map.
+    """
+    def __init__(self, *, ra1, ra2, dec1, dec2, value, nside_render=None):
+        self._ra1 = ra1
+        self._ra2 = ra2
+        self._dec1 = dec1
+        self._dec2 = dec2
+        self._value = value
+        self._nside_render = nside_render
+
+        sc_ra1 = np.isscalar(self._ra1)
+        sc_ra2 = np.isscalar(self._ra2)
+        sc_dec1 = np.isscalar(self._dec1)
+        sc_dec2 = np.isscalar(self._dec2)
+
+        if not sc_ra1 or not sc_ra2 or not sc_dec1 or not sc_dec2:
+            raise ValueError("Box only accepts scalar inputs for ra1, ra2, dec1, and dec2")
+
+    @property
+    def ra1(self):
+        return self._ra1
+
+    @property
+    def ra2(self):
+        return self._ra2
+
+    @property
+    def dec1(self):
+        return self._dec1
+
+    @property
+    def dec2(self):
+        return self._dec2
+
+    def _render(self, *, nside_render, return_pixel_ranges):
+        return hpg.query_box(
+            nside_render,
+            self._ra1,
+            self._ra2,
+            self._dec1,
+            self._dec2,
             nest=True,
             inclusive=False,
-            return_pixel_ranges=True,
+            return_pixel_ranges=return_pixel_ranges,
         )
 
     def __repr__(self):
-        s = 'Ellipse(ra=%.16g, dec=%16g, semi_major=%16g, semi_minor=%16g, alpha=%16g, value=%s)'
-        return s % (self._ra, self._dec, self._semi_major, self._semi_minor, self._alpha, repr(self._value))
+        s = "Box(ra1=%.16g, ra2=%.16g, dec1=%.16g, dec2=%.16g, value=%s, nside_render=%s"
+        return s % (self._ra1, self._ra2, self._dec1, self._dec2, repr(self._value), repr(self._nside_render))