diff --git a/lib/cartopy/crs.py b/lib/cartopy/crs.py
index 360ba93312..a9142b659f 100644
--- a/lib/cartopy/crs.py
+++ b/lib/cartopy/crs.py
@@ -1904,6 +1904,78 @@ def y_limits(self):
                                           semiminor_axis=6371007.181))
 
 
+class rHEALPix(Projection):
+    """
+    The rHEALPix projection is an extension of the HEALPix to present squares
+    (rather than triangles) at the north and south poles.
+
+    Parameters
+    ----------
+    central_longitude
+    north_square: int
+        The position for the north pole square. Must be one of 0, 1, 2 or 3.
+    south_square: int
+        The position for the south pole square. Must be one of 0, 1, 2 or 3.
+
+    """
+    def __init__(self, central_longitude=0, north_square=0, south_square=0):
+        valid_square = [0, 1, 2, 3]
+        if north_square not in valid_square:
+            raise ValueError('north_square must be one of {}'.format(valid_square))
+        if south_square not in valid_square:
+            raise ValueError('south_square must be one of {}'.format(valid_square))
+
+        proj4_params = [('proj', 'rhealpix'),
+                        ('north_square', north_square),
+                        ('south_square', south_square),
+                        ('lon_0', central_longitude)]
+        super(rHEALPix, self).__init__(proj4_params)
+
+        # Boundary is based on units of m, with a standard spherical ellipse.
+        # The hard-coded scale is the reason for not accepting the globe keyword.
+        nrth_x_pos = (north_square - 2) * 1e7
+        sth_x_pos = (south_square - 2) * 1e7
+        top = 5.05e6
+
+        points = [[2e7, -5e6],
+                  [2e7, top],
+                  [nrth_x_pos + 1e7, top],
+                  [nrth_x_pos + 1e7, 1.5e7],
+                  [nrth_x_pos, 1.5e7],
+                  [nrth_x_pos, top],
+                  [-2e7, top]]
+        if south_square != 0:
+            points.append([-2e7, -top])
+        points.extend([
+                  [sth_x_pos, -5e6],
+                  [sth_x_pos, -1.5e7],
+                  [sth_x_pos + 1e7, -1.5e7],
+                  [sth_x_pos + 1e7, -5e6],
+                  ])
+        self._boundary = sgeom.LineString(points[::-1])
+
+        xs, ys = zip(*points)
+        self._x_limits = min(xs), max(xs)
+        self._y_limits = min(ys), max(ys)
+        self._threshold = (self.x_limits[1] - self.x_limits[0]) / 1e4
+
+    @property
+    def boundary(self):
+        return self._boundary
+
+    @property
+    def threshold(self):
+        return self._threshold
+
+    @property
+    def x_limits(self):
+        return self._x_limits
+
+    @property
+    def y_limits(self):
+        return self._y_limits
+
+
 class _BoundaryPoint(object):
     def __init__(self, distance, kind, data):
         """
diff --git a/lib/cartopy/tests/test_crs.py b/lib/cartopy/tests/test_crs.py
index a9dc3ff413..acad1423f2 100644
--- a/lib/cartopy/tests/test_crs.py
+++ b/lib/cartopy/tests/test_crs.py
@@ -225,6 +225,24 @@ def test_utm(self):
                              decimal=1)
 
 
+def test_rHEALPix_defaults():
+    crs = ccrs.rHEALPix()
+    assert sorted(crs.proj4_params.items()) == [('ellps', 'WGS84'),
+                                                ('lon_0', 0),
+                                                ('north_square', 0),
+                                                ('proj', 'rhealpix'),
+                                                ('south_square', 0)]
+
+
+def test_rHEALPix_params():
+    crs = ccrs.rHEALPix(central_longitude=20, north_square=1, south_square=2)
+    assert sorted(crs.proj4_params.items()) == [('ellps', 'WGS84'),
+                                                ('lon_0', 20),
+                                                ('north_square', 1),
+                                                ('proj', 'rhealpix'),
+                                                ('south_square', 2)]
+
+
 def test_pickle():
     # check that we can pickle a simple CRS
     fh = BytesIO()