Extracting data from cartopy.feature

Question

how can I extract contour lines from data imported through cartopy's feature interface? If the solution involves geoviews.feature or another wrapper, that is OK, of course.

For instance, how would I extract the data plotted as cfeature.COASTLINE in the following example?

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

ax = plt.axes(projection=ccrs.PlateCarree())
ax.add_feature(cfeature.COASTLINE)
plt.show()

I'm grateful for any hints you might have!

FWIW, in basemap, I would do it like this:

import mpl_toolkits.basemap as bm
import matplotlib.pyplot as plt
m = bm.Basemap(width=2000e3,height=2000e3,
            resolution='l',projection='stere',
            lat_ts=70,lat_0=70,lon_0=-60.)

fig,ax=plt.subplots()
coastlines = m.drawcoastlines().get_segments()

Thomas Kühn · Accepted Answer · 2018-09-08T21:26:14.190

You can get the coordinates for the plotted lines directly from the feature, which contains a set of shapely.MultiLineStrings. As a proof of concept, check out this code:

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

fig, (ax1,ax2) = plt.subplots(nrows=2, subplot_kw = dict(projection=ccrs.PlateCarree()))
ax1.add_feature(cfeature.COASTLINE)

for geom in cfeature.COASTLINE.geometries():
    for g in geom.geoms:
        print(list(g.coords))
        ax2.plot(*zip(*list(g.coords)))

plt.show()

which gives this picture:

In other words, you can iterate over the MultiLineStrings of the feature by accessing its geometries(). Each of these MultiLineStrings then contains one or more LineStrings, which have a coords attribute that can be converted into a list. Hope this helps.

score 1 · Answer 2 · answered Sep 12 '18 at 10:51

For future reference: Some time later, I also came across this (more general?) method to access any feature:

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader

shpfilename = shpreader.natural_earth(resolution='110m',
                                      category='physical',
                                      name='coastline')  
coastlines = shpreader.Reader(shpfilename).records()

fig, ax = plt.subplots(subplot_kw = dict(projection=ccrs.PlateCarree()))
for c in coastlines:
    for g in c.geometry:
        ax.plot(*zip(*list(g.coords)))

yielding the same plot as above.

Extracting data from cartopy.feature

2 Answers2