如何使用xarray和cartopy在地图上绘制POP2网格（POP_gx1v7）数据？

Question

我正在处理 CESM2 的 LENS2 集合 的海洋变量，现在正在尝试在地理网格上绘制数据。海洋模型 POP2 基于极点位于格陵兰岛和南极洲的位移极点网格。 netCDF 文件包含相关的 WGS84 纬度和经度，但我无法绘制数据，特别是北纬 60 度以北的数据。 我的大部分工作都是使用

xarray

进行的，并且希望留在这个库中，但如果有必要，我也准备切换到其他工具。

理想情况下，我希望仍然能够看到原始的网格单元，而不是插入代码。 （如果我在北纬 60 度以下绘制数据，则可以使用

plt.contour()

进行插值。但是，我需要在北欧海域内绘制数据，并且该区域的图会崩溃。）

为了说明问题，这里是最简单的情况：

field = ds_temp.TEMP.isel(time=0, z_t=0)
extent = (ds_temp.TLONG.min(), ds_temp.TLONG.max(), 
          ds_temp.TLAT.min(), ds_temp.TLAT.max())

fig = plt.figure(figsize=(12, 6))
ax = plt.axes(projection=ccrs.PlateCarree())
cs = ax.imshow(field, extent=extent, 
               origin='lower', 
               transform=ccrs.PlateCarree())
ax.coastlines()
plt.show()

这会产生以下情节：

如您所见，

cartopy

海岸线和数据的陆地掩模发生了变化。我怎样才能用

data

绘制我的

cartopy

并如果可能的话继续使用

xarray

函数？

Answer 1

您可以尝试在大陆周围绘制自己的海岸线，在地图上显示为白色区域，而不是调用

ax.coastlines()

。最有可能的是，您在

field

中的原始数据对于地图的大陆区域有一些无效值。如果你打印它，它可能看起来像这样：

[[-1.7999999523162842 -1.7999999523162842 -1.7999999523162842 ...
  -1.7999999523162842 -1.7999999523162842 -1.7999999523162842]
 [-1.7999999523162842 -1.7999999523162842 -1.7999999523162842 ...
  -1.7999999523162842 -1.7999999523162842 -1.7999999523162842]
 [-1.7999999523162842 -1.7999999523162842 -1.7999999523162842 ...
  -1.7999999523162842 -1.7999999523162842 -1.7999999523162842]
 ...
 [-- -- -- ... -- -- --]
 [-- -- -- ... -- -- --]
 [-- -- -- ... -- -- --]]

那些双破折号看起来像是屏蔽数组中丢失的数据。因此，绘制海岸线的一种方法是使用 matplotlib.pyplot.contour 和稍微修改过的数据，其中这些“缺失值”被替换为其他一些值，允许您将它们与绘制轮廓的合适阈值进行比较：

import os
import matplotlib.pyplot as plt
from netCDF4 import Dataset as netcdf_dataset

from cartopy import config
import cartopy.crs as ccrs


# get the path of the file -i t can be found in the repo data directory.
fname = os.path.join(config["repo_data_dir"],
                     'netcdf', 'HadISST1_SST_update.nc'
                     )

# get the data, latitudes and longtitudes
dataset = netcdf_dataset(fname)
sst = dataset.variables['sst'][0, :, :]
lats = dataset.variables['lat'][:]
lons = dataset.variables['lon'][:]

ax = plt.axes(projection=ccrs.PlateCarree())

# plot data
ax.contourf(lons,
             lats,
             sst,
             60,
             transform=ccrs.PlateCarree()
             )

# replace invalid values (corresponding to land) with some unrealistically high values
very_high_value = 100
fixed_sst = sst.filled(very_high_value)

# use a slightly smaller threshold value to plot contour lines around land
ax.contour(lons,
           lats,
           fixed_sst,
           levels=[very_high_value-1],
           colors='black',
           transform=ccrs.PlateCarree()
           )

plt.show()

上面的代码基于这个例子，用matplotlib.pyplot.contourf展示数据，结果如下图：

看起来不是特别好看，但似乎可以实现你想要的功能。您需要选择合适的投影。如果您使用

imshow()

来绘制原始数据，那么您可能需要分别调整

extent

和

imshow()

中的

contour()

参数。