Creating a land/ocean mask in Cartopy?

Question

I have a gridded array of data with an associated lat/lon and I want to use Cartopy to find whether each particular grid cell is over ocean or land.

I can get the land data simply from the Cartopy Feature interface

land_10m = cartopy.feature.NaturalEarthFeature('physical', 'land', '10m')

I can then get the geometries

all_geometries = [geometry for geometry in land_50m.geometries()]

But I can't figure out how to use these geometries to test with a particular location is over land or not.

This question seems to have come up before Mask Ocean or Land from data using Cartopy and the solution was pretty much "use basemap instead" which is a bit unsatisfactory.

======== Update ========

Thanks to bjlittle I have a solution that works and generates the plot below

from IPython import embed
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
from shapely.geometry import Point
import cartopy

land_10m = cartopy.feature.NaturalEarthFeature('physical', 'land', '10m')
land_polygons = list(land_10m.geometries())

lats = np.arange(48,58, 0.1)
lons = np.arange(-5,5, 0.1)
lon_grid, lat_grid = np.meshgrid(lons, lats)

points = [Point(point) for point in zip(lon_grid.ravel(), lat_grid.ravel())]

land = []
for land_polygon in land_polygons:
    land.extend([tuple(point.coords)[0] for point in filter(land_polygon.covers, points)])

fig = plt.figure(figsize=(8, 8))
ax = fig.add_axes([0, 0, 1, 1], projection=ccrs.PlateCarree())

ax.add_feature(land_10m, zorder=0, edgecolor='black',     facecolor=sns.xkcd_rgb['black'])

xs, ys = zip(*land)
ax.scatter(xs, ys, transform=ccrs.PlateCarree(),
       s=12, marker='s', c='red', alpha=0.5, zorder=2)

plt.show()

However this is extremely slow and ultimately I'll need to do this globally at a much higher resolution.

Does anyone have any suggestions on how to adapt the above to be faster?

==== Update 2 ====

Preparing the polygons makes a HUGE difference

Adding these two lines has sped up an example at 1 degree from 40 seconds to 0.3 seconds

from shapely.prepared import prep
land_polygons_prep = [prep(land_polygon) for land_polygon in land_polygons]

The example that I ran at 0.1 degree that took over half an hour (i.e. it ran over lunch) now runs in 1.3 seconds :-)

Answer 1

Here's a code example using the solutions above that solves this problem:

from IPython import embed
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
from shapely.geometry import Point
import cartopy
from shapely.prepared import prep
import seaborn as sns


land_10m = cartopy.feature.NaturalEarthFeature('physical', 'land', '10m')
land_polygons = list(land_10m.geometries())

lats = np.arange(48,58, 0.1)
lons = np.arange(-5,5, 0.1)
lon_grid, lat_grid = np.meshgrid(lons, lats)

points = [Point(point) for point in zip(lon_grid.ravel(), lat_grid.ravel())]


land_polygons_prep = [prep(land_polygon) for land_polygon in land_polygons]


land = []
for land_polygon in land_polygons_prep:
    land.extend([tuple(point.coords)[0] for point in filter(land_polygon.covers, points)])

fig = plt.figure(figsize=(8, 8))
ax = fig.add_axes([0, 0, 1, 1], projection=ccrs.PlateCarree())

ax.add_feature(land_10m, zorder=0, edgecolor='black', facecolor=sns.xkcd_rgb['black'])

xs, ys = zip(*land)
ax.scatter(xs, ys, transform=ccrs.PlateCarree(),
       s=12, marker='s', c='red', alpha=0.5, zorder=2)

Answer 2

To highlight one approach that you could take, I've based my answer on the excellent Cartopy Hurricane Katrina (2005) gallery example, which plots the storm track of Katrina as it sweeps across the Gulf of Mexico, and up over the United States.

By essentially using a simple mixture of the cartopy.io.shapereader and some shapely predicates, we can get the job done. My example is a little bit verbose, but don't be put off... it's not that scary:

import matplotlib.pyplot as plt
import shapely.geometry as sgeom

import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader


def sample_data():
    """
    Return a list of latitudes and a list of longitudes (lons, lats)
    for Hurricane Katrina (2005).

    The data was originally sourced from the HURDAT2 dataset from AOML/NOAA:
    http://www.aoml.noaa.gov/hrd/hurdat/newhurdat-all.html on 14th Dec 2012.

    """
    lons = [-75.1, -75.7, -76.2, -76.5, -76.9, -77.7, -78.4, -79.0,
            -79.6, -80.1, -80.3, -81.3, -82.0, -82.6, -83.3, -84.0,
            -84.7, -85.3, -85.9, -86.7, -87.7, -88.6, -89.2, -89.6,
            -89.6, -89.6, -89.6, -89.6, -89.1, -88.6, -88.0, -87.0,
            -85.3, -82.9]

    lats = [23.1, 23.4, 23.8, 24.5, 25.4, 26.0, 26.1, 26.2, 26.2, 26.0,
            25.9, 25.4, 25.1, 24.9, 24.6, 24.4, 24.4, 24.5, 24.8, 25.2,
            25.7, 26.3, 27.2, 28.2, 29.3, 29.5, 30.2, 31.1, 32.6, 34.1,
            35.6, 37.0, 38.6, 40.1]

    return lons, lats


# create the figure and geoaxes
fig = plt.figure(figsize=(14, 12))
ax = fig.add_axes([0, 0, 1, 1], projection=ccrs.LambertConformal())
ax.set_extent([-125, -66.5, 20, 50], ccrs.Geodetic())

# load the shapefile geometries
shapename = 'admin_1_states_provinces_lakes_shp'
states_shp = shpreader.natural_earth(resolution='110m',
                                     category='cultural', name=shapename)
states = list(shpreader.Reader(states_shp).geometries())

# get the storm track lons and lats
lons, lats = sample_data()

# to get the effect of having just the states without a map "background"
# turn off the outline and background patches
ax.background_patch.set_visible(False)
ax.outline_patch.set_visible(False)

# turn the lons and lats into a shapely LineString
track = sgeom.LineString(zip(lons, lats))

# turn the lons and lats into shapely Points
points = [sgeom.Point(point) for point in zip(lons, lats)]

# determine the storm track points that fall on land
land = []
for state in states:
    land.extend([tuple(point.coords)[0] for point in filter(state.covers, points)])

# determine the storm track points that fall on sea
sea = set([tuple(point.coords)[0] for point in points]) - set(land)

# plot the state polygons
facecolor = [0.9375, 0.9375, 0.859375]    
for state in states:
    ax.add_geometries([state], ccrs.PlateCarree(),
                      facecolor=facecolor, edgecolor='black', zorder=1)

# plot the storm track land points 
xs, ys = zip(*land)
ax.scatter(xs, ys, transform=ccrs.PlateCarree(),
           s=100, marker='s', c='green', alpha=0.5, zorder=2)

# plot the storm track sea points
xs, ys = zip(*sea)
ax.scatter(xs, ys, transform=ccrs.PlateCarree(),
           s=100, marker='x', c='blue', alpha=0.5, zorder=2)

# plot the storm track
ax.add_geometries([track], ccrs.PlateCarree(),
                  facecolor='none', edgecolor='k')

plt.show()

The above example should generate the following plot, which highlights how to discriminate between land and sea points using shapely geometries.

HTH