寻找 Cartopy 选项在 zorder 不起作用时在数据上绘制土地

Question

GHRSST-MUR 卫星产品产生的海面温度在实际覆盖陆地的栅格的某些单元中具有数据。这种情况尤其发生在岛屿上。例如，如果您使用 ncview、NOAA ERDDAP 或 Cartopy 从 netCDF 绘制包括美国长岛在内的区域的地图，则海表温度数据将绘制在长岛上空，但不会绘制在大陆上。将 zorder 调整得更高以在分层中最后绘制 cartopy 土地要素并不能解决问题。我怀疑这是因为岛屿上的栅格单元包含该卫星产品中的数据。我正在 cartopy 中寻找一个选项来指定绘制发生位置的土地，即使数据数组在这些纬度/经度位置具有数据变量的数值。

我尝试为 cartopy land 特征设置 zorder=4 ，为 SST 的 pcolormesh 图设置 zorder=1 。这并没有按预期工作，因为陆地特征使海洋和陆地变灰，遮盖了海洋上空的海表温度数据。

我使用 ncview 绘制了相同的 netcdf 文件，以确认在长岛上绘制了 SST。 。

数据来源https://coastwatch.pfeg.noaa.gov/erddap/griddap/jplMURST41.nc?analysisd_sst%5B(2023-11-26T09:00:00Z)%5D%5B(38.2):(42.0)% 5D%5B(-74.2):(-67.9)%5D&.draw=表面&.vars=经度%7Clatitude%7Canalysis_sst&.colorBar=%7C%7C%7C%7C%7C&.bgColor=0xffccccff.

使用 NOAA ERDDAP 创建的图。请注意长岛上空的地罩是如何失效的。 https://coastwatch.pfeg.noaa.gov/erddap/griddap/jplMURST41.largePng?analysis_sst%5B(2023-11-26T09:00:00Z)%5D%5B(38.2):(42.0)%5D%5B (-74.2):(-67.9)%5D&.draw=surface&.vars=longitude%7Clatitude%7Canalysis_sst&.colorBar=%7C%7C%7C%7C%7C&.bgColor=0xffccccff

Answer 1

您应该分享您迄今为止尝试过的最小可重现示例。您提到设置 zorder 不起作用，但这正是您应该定义绘图顺序的方式。

使用像这样的常规网格，您还可以使用

imshow

，它允许仅使用角坐标指定范围，例如：

from osgeo import gdal
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

fn = "jplMURSST41_2ee1_731e_31b4.nc"

ds = gdal.OpenEx(fn)
data = ds.ReadAsArray()
gt = ds.GetGeoTransform()
ds = None

ys, xs = data.shape
ulx, uly = gdal.ApplyGeoTransform(gt, 0, 0)
lrx, lry = gdal.ApplyGeoTransform(gt, xs, ys)

mpl_extent = (ulx, lrx, lry, uly)

fig, ax = plt.subplots(
    figsize=(8,6), facecolor="w", layout="constrained",
    subplot_kw=dict(projection=ccrs.PlateCarree()),
)

im = ax.imshow(data, extent=mpl_extent, transform=ccrs.PlateCarree(), vmin=8, vmax=24, cmap="inferno", zorder=0)
cb = fig.colorbar(im, ax=ax, shrink=.4, label="SST [°C]")

ax.add_feature(cfeature.LAND, facecolor='w', zorder=1)

尝试交换 zorder 值，看看会发生什么。

使用

pcolormesh

或

pcolorfast

等类似，但这需要指定所有坐标。

# ...
ulx, uly = gdal.ApplyGeoTransform(gt, 0.5, 0.5)
lrx, lry = gdal.ApplyGeoTransform(gt, xs-0.5, ys-0.5)

xx = np.linspace(ulx, lrx, xs)
yy = np.linspace(uly, lry, ys)

# ...
im = ax.pcolorfast(xx, yy, data, transform=ccrs.PlateCarree(), vmin=8, vmax=24, cmap="inferno", zorder=0)
# ...

结果应该基本相同。当您将地图投影更改为与数据不同的东西时，这可能更重要。