如何在 Cartopy 的高分辨率海岸线上应用地球特征和陆地/海洋掩模？

Question

我在 Cartopy 中使用 GSHHS 数据集的海岸线。这对海岸线具有高分辨率。但我不仅要绘制高分辨率的海岸线，还要为海洋应用蒙版。

import matplotlib.pyplot as plt
import cartopy

fig = plt.figure(figsize=(20,12))
ax = plt.axes(projection=cartopy.crs.PlateCarree())
coast = cartopy.feature.GSHHSFeature(scale="full")
ax.add_feature(coast, linewidth=2)
ax.add_feature(cartopy.feature.NaturalEarthFeature("physical", "land", "10m"))
ax.set_extent([-17, -16, 27.9, 28.7])

执行代码时图像存在差异，因为我猜 ax.add_feature(cartopy.feature.NaturalEarthFeature("physical", "land", "10m")) 使用的是“10m”分辨率，而 GSHHS 具有更高的分辨率。

如何使用更高分辨率的 GSHHS 进行掩码？谢谢。

Answer 1

在回答问题how to apply a mask to hide features in the main plot之前，我们需要先调查可用的掩码。

在我们的例子中，main plot 是 Natural_Earth 10m 分辨率的 Physical Land 特征，GSHHSFeature 的各种分辨率作为可用掩码。

下面的代码和输出图揭示了这一点。

# Code adapted from:-
# Src: https://ctroupin.github.io/posts/2019-09-02-fine-coast/
import matplotlib.pyplot as plt
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature

resolutions = {"c": "crude",
               "l": "low",
               "i": "intermediate",
               "h": "high",
               "f": "full"}
coordinates = (8.7, 8.81, 42.55, 42.60)
myproj = ccrs.PlateCarree()
fig = plt.figure(figsize=(8, 4))
for i, res in enumerate(resolutions):
    ax = plt.subplot(2, 3, i+1, projection=myproj)
    coast = cfeature.GSHHSFeature(scale=res)
    ax.add_feature(coast, facecolor="lightgray")
    ax.add_feature(cartopy.feature.NaturalEarthFeature("physical", "land", "10m"),
              ec="red", fc="yellow", lw=2, alpha=0.4)
    ax.set_xlim(coordinates[0], coordinates[1])
    ax.set_ylim(coordinates[2], coordinates[3])
    plt.title(resolutions[res])
plt.suptitle("GSHHS: gray Versus 10m_Physical_Land: yellow/red")
plt.show()

假设我们需要一个在这个缩放级别的图。很明显，来自 2 个数据源的轮廓不太符合观众的眼睛。我们可能会得出结论，没有一个可用的掩码适合目标图。

但是，如果地块范围更宽，或者地块的比例更小，再加上一些制图技术，例如使用较厚的海岸线，可能会得到可接受的地块。该过程是反复试验的方法。

编辑1

添加 (Global_land_mask) 后，可以绘制更多选择 比较。

from global_land_mask import globe
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import numpy as np

# Extent of map in degrees
minlon,maxlon,minlat,maxlat = (8.7, 8.81, 42.55, 42.60)

# Lat/lon points to get for `global_land_mask` uses
# Finer than 500x250 has no improvement 
lons = np.linspace(minlon,maxlon, 500)
lats = np.linspace(minlat,maxlat, 250)

# Make a grid 
lon_grid, lat_grid = np.meshgrid(lons,lats)

# Get whether the points are on land.
z = globe.is_land(lat_grid, lon_grid)

# GSHHS ... 
resolutions = {"c": "crude",
               "l": "low",
               "i": "intermediate",
               "h": "high",
               "f": "full"}

myproj = ccrs.PlateCarree()
fig = plt.figure(figsize=(8, 4))
for i, res in enumerate(resolutions):
    ax = plt.subplot(2, 3, i+1, projection=myproj)

    # GSHHSFeature
    coast = cfeature.GSHHSFeature(scale=res)
    ax.add_feature(coast, facecolor="brown", alpha=0.5)

    # 10m physical_land
    ax.add_feature(cfeature.NaturalEarthFeature("physical", "land", "10m"),
              ec="red", fc="yellow", lw=2, alpha=0.4)

    # Global_land_mask data is used to create fillcontour
    # The fillcontour with proper (colormap, zorder, alpha) can be used as land `mask`
    ax.contourf(lon_grid, lat_grid, z, cmap="Greys_r", alpha=0.4)

    ax.set_xlim(minlon, maxlon)
    ax.set_ylim(minlat, maxlat)
    plt.title(resolutions[res])

plt.suptitle("GSHHS:brown/black | 10m_Land:yellow/red | Global_land_mask:light_gray")
plt.show()

# The best resolutuion from `Global_land_mask` is plotted in `lightgray` covering the sea areas