python - 重新网格化常规 netcdf 数据

Question

我有一个包含全局海面温度的 netcdf 文件。使用 matplotlib 和 Basemap，我设法使用以下代码制作了这些数据的 map :

from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

filename = '/Users/Nick/Desktop/SST/SST.nc'
fh = Dataset(filename, mode='r')

lons = fh.variables['LON'][:]
lats = fh.variables['LAT'][:]
sst = fh.variables['SST'][:].squeeze()

fig = plt.figure()

m = Basemap(projection='merc', llcrnrlon=80.,llcrnrlat=-25.,urcrnrlon=150.,urcrnrlat=25.,lon_0=115., lat_0=0., resolution='l')

lon, lat = np.meshgrid(lons, lats)
xi, yi = m(lon, lat)

cs = m.pcolormesh(xi,yi,sst, vmin=18, vmax=32)

m.drawmapboundary(fill_color='0.3')
m.fillcontinents(color='0.3', lake_color='0.3')
cbar = m.colorbar(cs, location='bottom', pad="10%", ticks=[18., 20., 22., 24., 26., 28., 30., 32.])
cbar.set_label('January SST (' + u'\u00b0' + 'C)')
plt.savefig('SST.png', dpi=300)

问题是数据的分辨率非常高(9km 网格)，这使得生成的图像非常嘈杂。我想将数据放在较低分辨率的网格上(例如 1 度)，但我正在努力弄清楚如何做到这一点。我遵循了一个可行的解决方案，通过将下面的代码插入到我上面的示例中来尝试使用 matplotlib griddata 函数，但它导致了“ValueError:条件必须是一维数组”。

xi, yi = np.meshgrid(lons, lats)

X = np.arange(min(x), max(x), 1)
Y = np.arange(min(y), max(y), 1)

Xi, Yi = np.meshgrid(X, Y)

Z = griddata(xi, yi, z, Xi, Yi)

我是 Python 和 matplotlib 的初学者，所以我不确定我做错了什么(或者更好的方法)。任何建议表示赞赏!

Answer 1

如果您使用例如双线性插值，这将导致更平滑场。

NCAR 气候数据指南有一个很好的 introduction to regridding (一般的，不是特定于 Python 的)。

据我所知，可用于 Python 的最强大的重新网格化例程实现是 Earth System Modeling Framework (ESMF) Python interface (ESMPy) .如果这对您的应用程序来说有点复杂，您应该查看

1. EarthPy重新网格化教程(例如使用 Pyresample 、 cKDTree 或 Basemap )。
2. 将您的数据转化为 Iris立方体和使用 Iris' regridding functions .

也许从查看 EarthPy regridding tutorial using Basemap 开始，因为您已经在使用它了。

在您的示例中执行此操作的方法是

from mpl_toolkits import basemap
from netCDF4 import Dataset

filename = '/Users/Nick/Desktop/SST/SST.nc'
with Dataset(filename, mode='r') as fh:
   lons = fh.variables['LON'][:]
   lats = fh.variables['LAT'][:]
   sst = fh.variables['SST'][:].squeeze()

lons_sub, lats_sub = np.meshgrid(lons[::4], lats[::4])

sst_coarse = basemap.interp(sst, lons, lats, lons_sub, lats_sub, order=1)

这会对您的 SST 数据执行双线性插值 (order=1) 到子采样网格(每四个点)。之后你的情节会看起来更粗粒度。如果您不喜欢那样，请使用例如

插值回原始网格

sst_smooth = basemap.interp(sst_coarse, lons_sub[0,:], lats_sub[:,0], *np.meshgrid(lons, lats), order=1)