3D CartoPy similar to Matplotlib-Basemap

Question

I'm new to Python with a question about Cartopy being able to be used in a 3D plot. Below is an example using matplotlibBasemap.

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.basemap import Basemap

m = Basemap(projection='merc',
            llcrnrlat=52.0,urcrnrlat=58.0,
            llcrnrlon=19.0,urcrnrlon=40.0,
            rsphere=6371200.,resolution='h',area_thresh=10)

fig = plt.figure()
ax = Axes3D(fig)
ax.add_collection3d(m.drawcoastlines(linewidth=0.25))
ax.add_collection3d(m.drawcountries(linewidth=0.35))
ax.add_collection3d(m.drawrivers(color='blue'))

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Height')

fig.show()

This creates a map within a 3D axis so that you can plot objects over the surface. But with Cartopy returns a matplotlib.axes.GeoAxesSubplot. Not clear how to take this and add to a 3D figure/axis as above with matplotlib-basemap.

So, can someone give any pointers on how to do a similar 3D plot with Cartopy?

Answer 1

The basemap mpl3d is a pretty neat hack, but it hasn't been designed to function in the described way. As a result, you can't currently use the same technique for much other than simple coastlines. For example, filled continents just don't work AFAICT.

That said, a similar hack is available when using cartopy. Since we can access shapefile information generically, this solution should work for any poly-line shapefile such as coastlines.

The first step is to get hold of the shapefile, and the respective geometries:

feature = cartopy.feature.NaturalEarthFeature('physical', 'coastline', '110m')
geoms = feature.geometries()

Next, we can convert these to the desired projection:

target_projection = ccrs.PlateCarree()
geoms = [target_projection.project_geometry(geom, feature.crs)
         for geom in geoms]

Since these are shapely geometries, we then want to convert them to matplotlib paths with:

from cartopy.mpl.patch import geos_to_path
import itertools

paths = list(itertools.chain.from_iterable(geos_to_path(geom)
                                             for geom in geoms))

With paths, we should be able to just create a PathCollection in matplotlib, and add it to the axes, but sadly, Axes3D doesn't seem to cope with PathCollection instances, so we need to workaround this by constructing a LineCollection (as basemap does). Sadly LineCollections don't take paths, but segments, which we can compute with:

segments = []
for path in paths:
    vertices = [vertex for vertex, _ in path.iter_segments()]
    vertices = np.asarray(vertices)
    segments.append(vertices)

Pulling this all together, we end up with a similar result to the basemap plot which your code produces:

import itertools

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection
import numpy as np

import cartopy.feature
from cartopy.mpl.patch import geos_to_path
import cartopy.crs as ccrs


fig = plt.figure()
ax = Axes3D(fig, xlim=[-180, 180], ylim=[-90, 90])
ax.set_zlim(bottom=0)


target_projection = ccrs.PlateCarree()

feature = cartopy.feature.NaturalEarthFeature('physical', 'coastline', '110m')
geoms = feature.geometries()

geoms = [target_projection.project_geometry(geom, feature.crs)
         for geom in geoms]

paths = list(itertools.chain.from_iterable(geos_to_path(geom) for geom in geoms))

# At this point, we start working around mpl3d's slightly broken interfaces.
# So we produce a LineCollection rather than a PathCollection.
segments = []
for path in paths:
    vertices = [vertex for vertex, _ in path.iter_segments()]
    vertices = np.asarray(vertices)
    segments.append(vertices)

lc = LineCollection(segments, color='black')

ax.add_collection3d(lc)

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Height')

plt.show()

On top of this, mpl3d seems to handle PolyCollection well, which would be the route I would investigate for filled geometries, such as the land outline (as opposed to the coastline, which is strictly an outline).

The important step is to convert the paths to polygons, and use these in a PolyCollection object:

concat = lambda iterable: list(itertools.chain.from_iterable(iterable))

polys = concat(path.to_polygons() for path in paths)
lc = PolyCollection(polys, edgecolor='black',
                    facecolor='green', closed=False)

The complete code for this case would look something like:

import itertools

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection, PolyCollection
import numpy as np

import cartopy.feature
from cartopy.mpl.patch import geos_to_path
import cartopy.crs as ccrs


fig = plt.figure()
ax = Axes3D(fig, xlim=[-180, 180], ylim=[-90, 90])
ax.set_zlim(bottom=0)


concat = lambda iterable: list(itertools.chain.from_iterable(iterable))

target_projection = ccrs.PlateCarree()

feature = cartopy.feature.NaturalEarthFeature('physical', 'land', '110m')
geoms = feature.geometries()

geoms = [target_projection.project_geometry(geom, feature.crs)
         for geom in geoms]

paths = concat(geos_to_path(geom) for geom in geoms)

polys = concat(path.to_polygons() for path in paths)

lc = PolyCollection(polys, edgecolor='black',
                    facecolor='green', closed=False)

ax.add_collection3d(lc)

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Height')

plt.show()

To yield: