ENH: Add HEALPix and rHEALPix projections

Author: greglucas

docs/make_projection.py

@@ -95,7 +95,8 @@ def utm_plot():
                   'OSGB': 4, 'LambertZoneII': 4.1, 'EuroPP': 5,
                   'Geostationary': 6, 'NearsidePerspective': 7,
                   'EckertI': 8.1, 'EckertII': 8.2, 'EckertIII': 8.3,
-                  'EckertIV': 8.4, 'EckertV': 8.5, 'EckertVI': 8.6}
+                  'EckertIV': 8.4, 'EckertV': 8.5, 'EckertVI': 8.6,
+                  'HEALPix': 9.1, 'RHEALPix': 9.2}
 
 
 def find_projections():

lib/cartopy/crs.py

@@ -1832,6 +1832,115 @@ def __init__(self):
         self.bounds = [-5242.32, 1212512.16, 1589155.51, 2706796.21]
 
 
+class HEALPix(Projection):
+    """
+    Hierarchical Equal Area isoLatitude Pixelisation (HEALPix) of a 2-sphere is an
+    algorithm for pixelisation of the 2-sphere based on subdivision of a distorted
+    rhombic dodecahedron. The HEALPix projection is area preserving and can be used
+    with a spherical and ellipsoidal model. It was initially developed for mapping
+    cosmic background microwave radiation.
+
+    """
+    _wrappable = True
+
+    def __init__(self, central_longitude=0, globe=None):
+        """
+        Parameters
+        ----------
+        central_longitude: optional
+            The true longitude of the central meridian in degrees.
+            Defaults to 0.
+
+        globe: optional
+            An instance of :class:`cartopy.crs.Globe`. If omitted, a default
+            globe with a spherical radius of 6370997 meters is created.
+        """
+        proj4_params = [('proj', 'healpix'),
+                        ('lon_0', central_longitude)]
+        super().__init__(proj4_params, globe=globe)
+        # Defaults to the PROJ sphere with semimajor axis 6370997m
+        w = (self.globe.semimajor_axis or 6370997) * np.pi
+        h = w/2
+        self.bounds = [-w, w, -h, h]
+        self._threshold = w / 1e4
+
+        self._boundary = sgeom.LinearRing(
+            [(-w, h/2), (-3/4*w, h), (-w/2, h/2), (-w/4, h), (0, h/2),
+             (w/4, h), (w/2, h/2), (3/4*w, h), (w, h/2),
+             (w, -h/2), (3/4*w, -h), (w/2, -h/2), (w/4, -h), (0, -h/2),
+             (-w/4, -h), (-w/2, -h/2), (-3/4*w, -h), (-w, -h/2), (-w, h/2)])
+
+    @property
+    def boundary(self):
+        return self._boundary
+
+
+class RHEALPix(Projection):
+    """
+    Also known as rHEALPix in PROJ, this projection is an extension of the
+    HEALPix projection to present rectangles, rather than triangles, at the
+    north and south poles.
+
+    Parameters
+    ----------
+    central_longitude: optional
+        The true longitude of the central meridian in degrees.
+        Defaults to 0.
+    north_square : int
+        The position for the north pole square. Must be one of 0, 1, 2 or 3.
+        0 would have the north pole square aligned with the left-most square,
+        and 3 would be aligned with the right-most.
+    south_square : int
+        The position for the south pole square. Must be one of 0, 1, 2 or 3.
+    """
+    _wrappable = True
+
+    def __init__(self, central_longitude=0, north_square=0, south_square=0, globe=None):
+        valid_square_pos = [0, 1, 2, 3]
+        if north_square not in valid_square_pos:
+            raise ValueError(f'north_square must be one of {valid_square_pos}')
+        if south_square not in valid_square_pos:
+            raise ValueError(f'south_square must be one of {valid_square_pos}')
+
+        proj4_params = [('proj', 'rhealpix'),
+                        ('north_square', north_square),
+                        ('south_square', south_square),
+                        ('lon_0', central_longitude)]
+        super().__init__(proj4_params)
+
+        # Defaults to the PROJ sphere with semimajor axis 6370997m
+        w = (self.globe.semimajor_axis or 6370997) * np.pi
+        h = 3/4 * w
+        box_h = w / 4
+        self.bounds = [-w, w, -h, h]
+        self._threshold = w / 1e4
+
+        self._boundary = sgeom.LinearRing([
+            # Left edge
+            (-w, box_h),
+            # Cut-out for the north square
+            (-w + north_square * w/2, box_h),
+            (-w + north_square * w/2, h),
+            (-w + (north_square + 1) * w/2, h),
+            (-w + (north_square + 1) * w/2, box_h),
+            # Right edge
+            (w, box_h),
+            (w, -box_h),
+            # Cut-out for the south square
+            (-w + (south_square + 1) * w/2, -box_h),
+            (-w + (south_square + 1) * w/2, -h),
+            (-w + south_square * w/2, -h),
+            (-w + south_square * w/2, -box_h),
+            # Left edge
+            (-w, -box_h),
+            (-w, box_h)
+            ])
+
+    @property
+    def boundary(self):
+        return self._boundary
+
+
 class LambertAzimuthalEqualArea(Projection):
     """
     A Lambert Azimuthal Equal-Area projection.

lib/cartopy/tests/crs/test_healpix.py

@@ -0,0 +1,22 @@
+# Copyright Crown and Cartopy Contributors
+#
+# This file is part of Cartopy and is released under the BSD 3-clause license.
+# See LICENSE in the root of the repository for full licensing details.
+"""
+Tests for the HEALPix projection.
+
+"""
+import cartopy.crs as ccrs
+from .helpers import check_proj_params
+
+
+def test_defaults():
+    crs = ccrs.HEALPix()
+    expected = {'ellps=WGS84', 'lon_0=0'}
+    check_proj_params('healpix', crs, expected)
+
+
+def test_central_longitude():
+    crs = ccrs.HEALPix(central_longitude=124.8)
+    expected = {'ellps=WGS84', 'lon_0=124.8'}
+    check_proj_params('healpix', crs, expected)

lib/cartopy/tests/crs/test_rhealpix.py

@@ -0,0 +1,36 @@
+# Copyright Crown and Cartopy Contributors
+#
+# This file is part of Cartopy and is released under the BSD 3-clause license.
+# See LICENSE in the root of the repository for full licensing details.
+"""
+Tests for the RHEALPix projection.
+
+"""
+import pytest
+
+import cartopy.crs as ccrs
+from .helpers import check_proj_params
+
+
+def test_defaults():
+    crs = ccrs.RHEALPix()
+    expected = {'ellps=WGS84', 'lon_0=0', 'north_square=0', 'south_square=0'}
+    check_proj_params('rhealpix', crs, expected)
+
+
+def test_square_positions():
+    crs = ccrs.RHEALPix(north_square=1, south_square=2)
+    expected = {'ellps=WGS84', 'lon_0=0', 'north_square=1', 'south_square=2'}
+    check_proj_params('rhealpix', crs, expected)
+
+def test_invalid_square_positions():
+    with pytest.raises(ValueError, match='north_square must be'):
+        ccrs.RHEALPix(north_square=4)
+    with pytest.raises(ValueError, match='south_square must be'):
+        ccrs.RHEALPix(south_square=-1)
+
+def test_central_longitude():
+    crs = ccrs.RHEALPix(north_square=2, central_longitude=-124.8)
+    expected = {'ellps=WGS84', 'lon_0=-124.8', 'north_square=2',
+                'south_square=0'}
+    check_proj_params('rhealpix', crs, expected)

lib/cartopy/tests/mpl/baseline_images/mpl/test_mpl_integration/test_global_map_HEALPix.png

@@ -169,6 +169,7 @@ def test_simple_global():
     ccrs.EqualEarth,
     ccrs.Gnomonic,
     ccrs.Hammer,
+    ccrs.HEALPix,
     pytest.param((ccrs.InterruptedGoodeHomolosine, dict(emphasis='land')),
                  id='InterruptedGoodeHomolosine'),
     ccrs.LambertCylindrical,
@@ -185,6 +186,7 @@ def test_simple_global():
     pytest.param((ccrs.RotatedPole,
                   dict(pole_latitude=45, pole_longitude=180)),
                  id='RotatedPole'),
+    ccrs.RHEALPix,
     ccrs.Stereographic,
     ccrs.SouthPolarStereo,
     pytest.param((ccrs.TransverseMercator, dict(approx=True)),