Projection extension - Gnomonic cubed-sphere

[darothen]

I have some upcoming work where I intend to develop some visualization tools for model output on a cubed-sphere grid, and I'd really like to use cartopy.

To gear up in this direction, I've been playing with some sample output data from NASA's GEOS-5 to get a sense of what grid specification and metadata might be useful to have and where friction points might be. Plotting individual faces of the output data work relatively well (although I do run into issues with pcolormesh wrapping around the international dateline for one of the faces). But ideally, I need to visualize my entire output on a single projection. The orthographic projection works ok, but only for a subset of the globe. Instead, a better way would be to "unfold" the cubed sphere and arrange the tiles in particular patterns.

How would one go about building such an "unfolded" projection into cartopy? I've toyed around with creating a grid of axes and plotting whatever faces may partially overlap on them, and that works okay, but it's extremely difficult to line up longitude lines with the "branched" tiles. For an example of how this could work, here's a sample plot made with Panoply:

Anyways. I'm hoping to do most of the implementation and coding myself. But if there are any resources/notes available for building new projections and extending proj.4, that would be extremely helpful!

[QuLogic]

The first place I'd go look is Proj.4 to see if they have the projection. I don't see anything that looks like a "Gnomonic cubed-sphere" at first glance, but there are some similar ones like rHEALPix and Quadrilateralized Spherical Cube.

Neither of these (nor any you might add to proj4, naturally) are available in Cartopy at the moment, but Proj.4 projections are rather straightforward to add. You can look at this Azimuthal Equidistant PR or this Alber's Equal Area PR for inspiration there.

[pelson]

I decided to have a go at implementing the rHEALPix projection as it has an inverse transformation defined (a pre-requisite of a cartopy.crs.Projection to be fully functional). Once I'd figured out the magic numbers for the boundaries, it turned out to be fairly painless:

from __future__ import division
from cartopy.crs import Projection
import shapely.geometry as sgeom

import cartopy.crs as ccrs


class rHEALPix(Projection):
    def __init__(self, central_longitude=0, north_square=0, south_square=0):
        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.
        nrth_x_pos = (north_square - 2) * 1e7
        sth_x_pos = (south_square - 2) * 1e7
        top = 5.05e6
        points = []
        points.extend([
                  [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

Naturally, there is a bunch of input validation + code hygeine that would need to be taken care of for this to become part of cartopy, but with the above (and a separate bugfix to cartopy that I'm about to make a PR for) the following code functions nicely:

import matplotlib.pyplot as plt
ax = plt.axes(projection=rHEALPix())
ax.stock_img()
ax.coastlines()
ax.gridlines()
ax.set_global()
plt.show()

and

import cartopy.crs as ccrs
import cartopy.feature
import matplotlib.pyplot as plt

ax = plt.axes(projection=rHEALPix(north_square=1, south_square=3, central_longitude=180))
ax.coastlines()
ax.gridlines()
ax.add_feature(cartopy.feature.LAND)
ax.add_feature(cartopy.feature.OCEAN)
plt.plot([50, -0.08, 132], [-30, 51.53, 43.17], color='red',
         transform=ccrs.Geodetic(), linewidth=3)
ax.set_global()
plt.show()

---

Note: The bugfix was relating to the use of error code 15 in proj:

diff --git a/lib/cartopy/_crs.pyx b/lib/cartopy/_crs.pyx
index 95b18ae..6795ead 100644
--- a/lib/cartopy/_crs.pyx
+++ b/lib/cartopy/_crs.pyx
@@ -282,7 +282,7 @@ cdef class CRS:
             cy *= DEG_TO_RAD
         status = pj_transform(src_crs.proj4, self.proj4, 1, 1, &cx, &cy, NULL);
 
-        if trap and status == -14 or status == -20:
+        if trap and status in [-14, -15, -20]:
             # -14 => "latitude or longitude exceeded limits"
             # -20 => "tolerance condition error"
             cx = cy = NAN
diff --git a/lib/cartopy/_trace.cpp b/lib/cartopy/_trace.cpp
index c8d2f4c..c5ef135 100644
--- a/lib/cartopy/_trace.cpp
+++ b/lib/cartopy/_trace.cpp
@@ -63,7 +63,7 @@ Point CartesianInterpolator::project(const Point &src_xy)
     xy.v = src_xy.y;
 
     int status = pj_transform(m_src_proj, m_dest_proj, 1, 1, &xy.u, &xy.v, NULL);
-    if (status == -14 || status == -20)
+    if (status == -14 || status == -20 || status == -15)
     {
         // -14 => "latitude or longitude exceeded limits"
         // -20 => "tolerance condition error"

[pelson]

A quick hack also produces the full healpix projection:

from __future__ import division

import matplotlib.pyplot as plt

import cartopy.crs as ccrs


from cartopy.crs import Projection
import shapely.geometry as sgeom

import cartopy.crs as ccrs


class Healpix(Projection):
    def __init__(self, central_longitude=0):
        proj4_params = [('proj', 'healpix'),
                        ('lon_0', central_longitude)]
        super(Healpix, self).__init__(proj4_params)

        # Boundary is based on units of m, with a standard spherical ellipse.
        width = 2e7
        h = width/2
        box_h = width/4
 
        points = [[width, -box_h],
                  [width, box_h],
                  [width - 1*h + h/2, h],
                  [width - 1*h, box_h],
                  [width - 2*h + h / 2, h],
                  [width - 2*h, box_h],
                  [width - 3*h + h / 2, h],
                  [width - 3*h, box_h],
                  [width - 4*h + h / 2, h],
                  [width - 4*h, box_h],
                  [width - 4*h, -box_h],
                  [width - 4*h + h / 2, -h],
                  [width - 3*h, -box_h],
                  [width - 3*h + h / 2, -h],
                  [width - 2*h, -box_h],
                  [width - 2*h + h / 2, -h],
                  [width - 1*h, -box_h],
                  [width - 1*h + h / 2, -h]]
        self._boundary = sgeom.LinearRing(points)

        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


ax = plt.axes(projection=Healpix(central_longitude=30))
ax.stock_img()
ax.coastlines()
plt.show()

Unfortunately it looks like there is a bug with proj4 (checked with 4.8.0 through to master) which puts the north polar triangles in the wrong place...