ENH: Add HEALPix and rHEALPix projections

Author: greglucas

docs/make_projection.py

@@ -6,9 +6,11 @@
 import inspect
 from pathlib import Path
 import textwrap
+import warnings
 
 import matplotlib.pyplot as plt
 import numpy as np
+from pyproj.exceptions import CRSError
 
 import cartopy.crs as ccrs
 
@@ -96,7 +98,8 @@ def utm_plot():
                   'Geostationary': 6, 'NearsidePerspective': 7,
                   'EckertI': 8.1, 'EckertII': 8.2, 'EckertIII': 8.3,
                   'EckertIV': 8.4, 'EckertV': 8.5, 'EckertVI': 8.6,
-                  'Spilhaus': 9}
+                  'Spilhaus': 9,
+                  'HEALPix': 9.1, 'RHEALPix': 9.2}
 
 
 def find_projections():
@@ -162,7 +165,14 @@ def prj_class_sorter(cls):
             # Get instance arguments and number of plots
             instance_args = SPECIFIC_PROJECTION_KWARGS.get(prj, {})
 
-            prj_inst, instance_repr = create_instance(prj, instance_args)
+            try:
+                prj_inst, instance_repr = create_instance(prj, instance_args)
+            except CRSError:
+                # Installed PROJ may be too old/new to support this projection
+                table.write('\n')
+                warnings.warn(
+                    f'Projection {prj.__name__} not supported by this version of PROJ')
+                continue
 
             aspect = (np.diff(prj_inst.x_limits) /
                       np.diff(prj_inst.y_limits))[0]

docs/source/reference/projections.rst

@@ -509,13 +509,34 @@ Spilhaus
 
 .. autoclass:: cartopy.crs.Spilhaus
 
+HEALPix
+-------
+
+.. autoclass:: cartopy.crs.HEALPix
+
+.. plot::
+
+    import matplotlib.pyplot as plt
+    import cartopy.crs as ccrs
+
+    plt.figure(figsize=(6, 3))
+    ax = plt.axes(projection=ccrs.HEALPix())
+    ax.coastlines(resolution='110m')
+    ax.gridlines()
+
+
+RHEALPix
+--------
+
+.. autoclass:: cartopy.crs.RHEALPix
+
 .. plot::
 
     import matplotlib.pyplot as plt
     import cartopy.crs as ccrs
 
-    plt.figure(figsize=(3.0004, 3))
-    ax = plt.axes(projection=ccrs.Spilhaus())
+    plt.figure(figsize=(4, 3))
+    ax = plt.axes(projection=ccrs.RHEALPix())
     ax.coastlines(resolution='110m')
     ax.gridlines()
 

lib/cartopy/crs.py

@@ -1862,6 +1862,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

@@ -173,6 +173,7 @@ def test_simple_global():
     ccrs.EqualEarth,
     ccrs.Gnomonic,
     ccrs.Hammer,
+    ccrs.HEALPix,
     pytest.param((ccrs.InterruptedGoodeHomolosine, dict(emphasis='land')),
                  id='InterruptedGoodeHomolosine'),
     ccrs.LambertCylindrical,
@@ -193,6 +194,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)),