Создание маски земли / океана в Cartopy?

Question

У меня есть массив данных с привязкой к сетке со связанным широтой / долготой, и я хочу использовать Cartopy, чтобы определить, находится ли каждая конкретная ячейка сетки над океаном или над сушей.

Я могу получить данные о земле просто изИнтерфейс Cartopy Feature

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

Затем я могу получить геометрию

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

Но я не могу понять, как использовать эти геометрии для тестирования в определенном месте над землей или нет.

Похоже, этот вопрос возник раньше, чем Маска Океана или Земли из данных с использованием 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

Вот пример кода, использующий решения, описанные выше, который решает эту проблему:

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

Чтобы выделить один подход, который вы могли бы предпринять, я основываю свой ответ на превосходном примере галереи Cartopy Hurricane Katrina (2005) , который показывает штормовой след Катрины какон охватывает Мексиканский залив и США.

По сути, используя простую смесь cartopy.io.shapereader и некоторых фигурных предикатов мы можем сделать работу.Мой пример немного многословен, но не стоит откладывать ... это не так страшно:

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

Приведенный выше пример должен сгенерировать следующий сюжет , который подчеркивает, какразличать сухопутные и морские точки, используя стройную геометрию.

HTH