diff --git a/docs/make_projection.py b/docs/make_projection.py index 8f50fd559..b0aca511c 100644 --- a/docs/make_projection.py +++ b/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(): diff --git a/lib/cartopy/crs.py b/lib/cartopy/crs.py index 600f8ac24..d3e4fbd7b 100644 --- a/lib/cartopy/crs.py +++ b/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. diff --git a/lib/cartopy/tests/crs/test_healpix.py b/lib/cartopy/tests/crs/test_healpix.py new file mode 100644 index 000000000..01c1ebb0d --- /dev/null +++ b/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) diff --git a/lib/cartopy/tests/crs/test_rhealpix.py b/lib/cartopy/tests/crs/test_rhealpix.py new file mode 100644 index 000000000..552488ddc --- /dev/null +++ b/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) diff --git a/lib/cartopy/tests/mpl/baseline_images/mpl/test_mpl_integration/test_global_map_HEALPix.png b/lib/cartopy/tests/mpl/baseline_images/mpl/test_mpl_integration/test_global_map_HEALPix.png new file mode 100644 index 000000000..273f4d6ab Binary files /dev/null and b/lib/cartopy/tests/mpl/baseline_images/mpl/test_mpl_integration/test_global_map_HEALPix.png differ diff --git a/lib/cartopy/tests/mpl/baseline_images/mpl/test_mpl_integration/test_global_map_RHEALPix.png b/lib/cartopy/tests/mpl/baseline_images/mpl/test_mpl_integration/test_global_map_RHEALPix.png new file mode 100644 index 000000000..17c836c76 Binary files /dev/null and b/lib/cartopy/tests/mpl/baseline_images/mpl/test_mpl_integration/test_global_map_RHEALPix.png differ diff --git a/lib/cartopy/tests/mpl/test_mpl_integration.py b/lib/cartopy/tests/mpl/test_mpl_integration.py index 6b528b982..6b0f0a9b9 100644 --- a/lib/cartopy/tests/mpl/test_mpl_integration.py +++ b/lib/cartopy/tests/mpl/test_mpl_integration.py @@ -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)),