python - 在Cartopy中创建陆地/海洋 mask ？

我有一个带有相关纬度/经度的网格数据数组，我想使用Cartopy来查找每个特定的网格单元是在海洋还是陆地上。
我可以简单地从cartpy特性界面获得土地数据

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

我可以得到几何图形

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

但我不知道如何使用这些几何图形来测试特定位置是否在陆地上。
这个问题似乎在Mask Ocean or Land from data using Cartopy之前就出现了，解决方案几乎是“改用basemap”，这有点不令人满意。
====更新========
多亏了bjlittle，我有了一个可行的解决方案，并生成了下面的图

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

但是这是非常缓慢的，最终我需要在全球范围内以更高的分辨率来完成。
有没有人对如何使上述方法更快速有什么建议？
===更新2====
准备多边形有很大的不同
把这两条线相加，使一个1度的例子从40秒加速到0.3秒

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

我在0.1度的温度下跑了半个多小时（也就是在午餐时间跑了），现在只需1.3秒

最佳答案

为了强调一种你可以采用的方法，我把我的答案放在了优秀的CARPOTY Hurricane Katrina (2005)画廊的例子上，它描绘了卡特丽娜横扫墨西哥湾和美国上空的风暴轨迹。
本质上，通过使用cartopy.io.shapereader和一些shapely predicates的简单混合，我们可以完成这项工作。我的例子有点冗长，但不要拖延。。。没那么可怕：

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

上面的例子应该生成以下plot，它强调了如何使用形状几何来区分陆地和海洋点。
高温高压

关于python - 在Cartopy中创建陆地/海洋 mask ？，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/53322952/