在 Cartopy 中创建陆地/海洋面具？

Question

我有一个带有关联纬度/经度的网格数据数组，我想使用 Cartopy 来查找每个特定的网格单元是在海洋上还是在陆地上。

我可以简单地从 Cartopy Feature 界面获取土地数据

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 之前提出的，解决方案几乎是“改用底图”，这有点不令人满意。

======== 更新========

感谢 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 秒 :-)

Answer 1

为了突出您可以采用的一种方法，我的回答基于出色的 Cartopy 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，它强调了如何使用有形状的几何图形来区分陆地和海洋点。

HTH

Answer 2

这是一个使用上述解决方案解决此问题的代码示例：

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)