Vfrdtyky

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

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

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)