diff --git a/docs/geometry.rst b/docs/geometry.rst index 64fac92..e596942 100644 --- a/docs/geometry.rst +++ b/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 diff --git a/healsparse/__init__.py b/healsparse/__init__.py index d9a0233..44f0544 100644 --- a/healsparse/__init__.py +++ b/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 diff --git a/healsparse/geom.py b/healsparse/geom.py index 2c113e3..566e9b9 100644 --- a/healsparse/geom.py +++ b/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,9 +447,9 @@ 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, @@ -417,21 +457,85 @@ def get_pixels(self, *, nside): 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)) diff --git a/tests/test_geom.py b/tests/test_geom.py index 3e028b2..de1018a 100644 --- a/tests/test_geom.py +++ b/tests/test_geom.py @@ -2,9 +2,10 @@ import numpy.testing as testing import numpy as np +import hpgeom as hpg import healsparse -from healsparse import Circle, Polygon, Ellipse +from healsparse import Circle, Polygon, Ellipse, Box def atbound(longitude, minval, maxval): @@ -201,6 +202,45 @@ def test_bad_values(self): value=1.0, ) + def test_circle_nside_render(self): + """Test using a circle with a different rendering nside.""" + nside = 2**17 + + ra, dec = 200.0, 0.0 + radius = 30.0/3600.0 + nside_render = 2**14 + circle = Circle( + ra=ra, + dec=dec, + radius=radius, + value=2**4, + nside_render=nside_render, + ) + + pixels = circle.get_pixels(nside=nside) + + pixels_coarse = hpg.query_circle(nside_render, ra, dec, radius, inclusive=False) + pixels_fine = hpg.upgrade_pixels(nside_render, pixels_coarse, nside) + + np.testing.assert_array_equal(pixels, pixels_fine) + + pixel_ranges = circle.get_pixel_ranges(nside=nside) + + pixel_ranges_coarse = hpg.query_circle( + nside_render, + ra, + dec, + radius, + inclusive=False, + return_pixel_ranges=True, + ) + pixel_ranges_fine = hpg.upgrade_pixel_ranges(nside_render, pixel_ranges_coarse, nside) + + np.testing.assert_array_equal(pixel_ranges, pixel_ranges_fine) + + with self.assertRaises(ValueError): + circle.get_pixels(nside=2**10) + def test_polygon_smoke(self): """ just test we can make a polygon and a map from it @@ -283,6 +323,57 @@ def test_polygon_values(self): testing.assert_array_equal(smap.valid_pixels, smap2.valid_pixels) testing.assert_array_equal(smap2.get_values_pix(smap2.valid_pixels), 2.0) + # Test booleans + poly = Polygon( + ra=ra, + dec=dec, + value=True, + ) + smap3 = poly.get_map(nside_coverage=32, nside_sparse=nside, dtype=np.bool_) + testing.assert_array_equal(smap.valid_pixels, smap3.valid_pixels) + testing.assert_array_equal(smap3.get_values_pix(smap3.valid_pixels), True) + + def test_polygon_nside_render(self): + """Test using a polygon with a different rendering nside.""" + nside = 2**17 + + # make a box + ra_range = 200.0, 200.1 + dec_range = 0.1, 0.2 + + ra = [ra_range[0], ra_range[1], ra_range[1], ra_range[0]] + dec = [dec_range[0], dec_range[0], dec_range[1], dec_range[1]] + nside_render = 2**14 + poly = Polygon( + ra=ra, + dec=dec, + value=64, + nside_render=nside_render, + ) + + pixels = poly.get_pixels(nside=nside) + + pixels_coarse = hpg.query_polygon(nside_render, ra, dec, inclusive=False) + pixels_fine = hpg.upgrade_pixels(nside_render, pixels_coarse, nside) + + np.testing.assert_array_equal(pixels, pixels_fine) + + pixel_ranges = poly.get_pixel_ranges(nside=nside) + + pixel_ranges_coarse = hpg.query_polygon( + nside_render, + ra, + dec, + inclusive=False, + return_pixel_ranges=True, + ) + pixel_ranges_fine = hpg.upgrade_pixel_ranges(nside_render, pixel_ranges_coarse, nside) + + np.testing.assert_array_equal(pixel_ranges, pixel_ranges_fine) + + with self.assertRaises(ValueError): + poly.get_pixels(nside=2**10) + def test_ellipse_smoke(self): """ Test that we can make an ellipse and a map from it. @@ -375,6 +466,133 @@ def test_ellipse_values(self): testing.assert_array_equal(smap.valid_pixels, smap2.valid_pixels) testing.assert_array_equal(smap2.get_values_pix(smap2.valid_pixels), 2.0) + def test_ellipse_render(self): + """Test using an ellipse with a different rendering nside.""" + nside = 2**17 + + ra, dec = 200.0, 0.0 + semi_major = 30.0/3600.0 + semi_minor = 15.0/3600.0 + alpha = 45.0 + nside_render = 2**14 + ellipse = Ellipse( + ra=ra, + dec=dec, + semi_major=semi_major, + semi_minor=semi_minor, + alpha=alpha, + value=2**4, + nside_render=nside_render, + ) + + pixels = ellipse.get_pixels(nside=nside) + + pixels_coarse = hpg.query_ellipse( + nside_render, + ra, + dec, + semi_major, + semi_minor, + alpha, + inclusive=False, + ) + pixels_fine = hpg.upgrade_pixels(nside_render, pixels_coarse, nside) + + np.testing.assert_array_equal(pixels, pixels_fine) + + pixel_ranges = ellipse.get_pixel_ranges(nside=nside) + + pixel_ranges_coarse = hpg.query_ellipse( + nside_render, + ra, + dec, + semi_major, + semi_minor, + alpha, + inclusive=False, + return_pixel_ranges=True, + ) + pixel_ranges_fine = hpg.upgrade_pixel_ranges(nside_render, pixel_ranges_coarse, nside) + + np.testing.assert_array_equal(pixel_ranges, pixel_ranges_fine) + + with self.assertRaises(ValueError): + ellipse.get_pixels(nside=2**10) + + def test_box_smoke(self): + """Test that we can make a box and a map from it.""" + ra1, ra2 = 20.0, 21.0 + dec1, dec2 = 5.0, 6.0 + nside = 2**17 + box = Box( + ra1=ra1, + ra2=ra2, + dec1=dec1, + dec2=dec2, + value=2**4, + ) + + pixels = box.get_pixels(nside=nside) + self.assertGreater(pixels.size, 0) + + smap = box.get_map(nside_coverage=32, nside_sparse=nside, dtype=np.int16) + self.assertTrue(isinstance(smap, healsparse.HealSparseMap)) + + smap2 = box.get_map_like(smap) + self.assertEqual(smap2.nside_coverage, smap.nside_coverage) + self.assertEqual(smap2.nside_sparse, smap.nside_sparse) + self.assertEqual(smap2.dtype, smap.dtype) + + def test_box_bad_values(self): + """Test that we can only make boxes with scalars.""" + ra1, ra2 = (20.0, 20.5), 21.0 + dec1, dec2 = 5.0, 6.0 + + with self.assertRaises(ValueError): + Box(ra1=ra1, ra2=ra2, dec1=dec1, dec2=dec2, value=0) + + def test_box_nside_render(self): + """Test using a box with a different rendering nside.""" + nside = 2**17 + + ra1, ra2 = 20.0, 21.0 + dec1, dec2 = 5.0, 6.0 + + nside_render = 2**14 + box = Box( + ra1=ra1, + ra2=ra2, + dec1=dec1, + dec2=dec2, + value=2**4, + nside_render=nside_render, + ) + + pixels = box.get_pixels(nside=nside) + + pixels_coarse = hpg.query_box(nside_render, ra1, ra2, dec1, dec2, inclusive=False) + pixels_fine = hpg.upgrade_pixels(nside_render, pixels_coarse, nside) + + np.testing.assert_array_equal(pixels, pixels_fine) + + pixel_ranges = box.get_pixel_ranges(nside=nside) + + pixel_ranges_coarse = hpg.query_box( + nside_render, + ra1, + ra2, + dec1, + dec2, + inclusive=False, + return_pixel_ranges=True, + ) + pixel_ranges_fine = hpg.upgrade_pixel_ranges(nside_render, pixel_ranges_coarse, nside) + + np.testing.assert_array_equal(pixel_ranges, pixel_ranges_fine) + + with self.assertRaises(ValueError): + box.get_pixels(nside=2**10) + def test_realize_geom_or(self): """ test "or"ing two geom objects