Search code examples
pythonmatplotlib-basemapcartopygeoviews

Extracting data from cartopy.feature


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()

Solution

  • 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:

    result of the above code

    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.