Création d'un masque terre/océan dans Cartopy ?

Question

J'ai un tableau de données en grille avec une latitude/longitude associée et je veux utiliser Cartopy pour déterminer si chaque cellule de grille particulière est sur l'océan ou la terre.

Je peux obtenir simplement les données terrestres depuis l'interface de Cartopy

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

Je peux ensuite obtenir les géométries

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

Mais je ne sais pas comment utiliser ces géométries pour tester si un emplacement particulier est sur la terre ou non.

Cette question semble avoir déjà été posée Mask Ocean or Land from data using Cartopy et la solution était en gros "utilisez basemap à la place", ce qui est un peu insatisfaisant.

\======== Mise à jour ========

Merci à bjlittle, j'ai une solution qui fonctionne et génère le graphique ci-dessous

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

Cependant, c'est extrêmement lent et finalement je devrai le faire à l'échelle mondiale à une résolution beaucoup plus élevée.

Est-ce que quelqu'un a des suggestions sur la façon d'adapter ce qui précède pour être plus rapide?

\==== Mise à jour 2 ====

Préparer les polygones fait une ÉNORME différence

Ajouter ces deux lignes a accéléré un exemple à 1 degré de 40 secondes à 0,3 seconde

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

L'exemple que j'ai exécuté à 0,1 degré et qui a pris plus d'une demi-heure (c'est-à-dire qu'il s'est exécuté pendant le déjeuner) se termine maintenant en 1,3 seconde :-)

Answer 1

Pour mettre en avant une approche que vous pourriez adopter, j'ai basé ma réponse sur l'excellent exemple de la galerie Cartopy Ouragan Katrina (2005), qui trace la trajectoire de l'ouragan Katrina alors qu'il balaye le golfe du Mexique et remonte sur les États-Unis.

En utilisant essentiellement un simple mélange du cartopy.io.shapereader et de certains prédispositions shapely, nous pouvons faire le travail. Mon exemple est un peu verbeux, mais ne soyez pas découragés... ce n'est pas si effrayant:

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

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

def sample_data():
    """
    Retourne une liste de latitudes et une liste de longitudes (lons, lats)
    pour l'ouragan Katrina (2005).

    Les données proviennent initialement de l'ensemble de données HURDAT2 de AOML/NOAA:
    http://www.aoml.noaa.gov/hrd/hurdat/newhurdat-all.html le 14 décembre 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

# créer la figure et les 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())

# charger les géométries du fichier shapefile
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())

# obtenir les longitudes et latitudes de la trajectoire de l'ouragan
lons, lats = sample_data()

# pour obtenir l'effet d'avoir seulement les états sans un arrière-plan de carte
# désactiver les contours et les fonds de carte
ax.background_patch.set_visible(False)
ax.outline_patch.set_visible(False)

# transformer les longitudes et latitudes en LineString shapely
track = sgeom.LineString(zip(lons, lats))

# transformer les longitudes et latitudes en Points shapely
points = [sgeom.Point(point) for point in zip(lons, lats)]

# déterminer les points de la trajectoire de l'ouragan qui se trouvent sur la terre
land = []
for state in states:
    land.extend([tuple(point.coords)[0] for point in filter(state.covers, points)])

# déterminer les points de la trajectoire de l'ouragan qui se trouvent en mer
sea = set([tuple(point.coords)[0] for point in points]) - set(land)

# tracer les polygones des états
facecolor = [0.9375, 0.9375, 0.859375]    
for state in states:
    ax.add_geometries([state], ccrs.PlateCarree(),
                      facecolor=facecolor, edgecolor='black', zorder=1)

# tracer les points de la trajectoire de l'ouragan sur terre
xs, ys = zip(*land)
ax.scatter(xs, ys, transform=ccrs.PlateCarree(),
           s=100, marker='s', c='green', alpha=0.5, zorder=2)

# tracer les points de la trajectoire de l'ouragan en mer
xs, ys = zip(*sea)
ax.scatter(xs, ys, transform=ccrs.PlateCarree(),
           s=100, marker='x', c='blue', alpha=0.5, zorder=2)

# tracer la trajectoire de l'ouragan
ax.add_geometries([track], ccrs.PlateCarree(),
                  facecolor='none', edgecolor='k')

plt.show()

L'exemple ci-dessus devrait générer le graphique suivant, qui illustre comment différencier les points de terre et de mer en utilisant les géométries shapely.

HTH

Answer 2

Voici un exemple de code utilisant les solutions ci-dessus qui résout ce problème :

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)