重网格化常规NetCDF数据
我有一个netcdf文件，里面包含了全球海面温度的数据。通过使用matplotlib和Basemap，我成功地用以下代码制作了这个数据的地图：
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)
问题是，这些数据的分辨率非常高（每个格子9公里），导致生成的图像看起来很杂乱。我想把数据放到一个更低分辨率的网格上（比如1度），但我不知道该怎么做。我尝试按照一个解决方案，使用matplotlib的griddata函数，把下面的代码插入到我之前的例子中，但结果出现了'ValueError: condition must be a 1-d array'的错误。
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还是个相对初学者，所以我不太确定自己哪里出错了（或者有什么更好的方法）。任何建议都非常感谢！
5 个回答
看看这个关于xarray的例子……使用ds.interp这个方法，并指定新的纬度和经度值。
http://xarray.pydata.org/en/stable/interpolation.html#example
如果你在使用Linux系统，可以通过nctoolkit这个工具来实现这个功能（https://nctoolkit.readthedocs.io/en/latest/）。
你没有说明你的数据的经纬度范围，所以我假设你是在处理全球的数据。将数据调整到1度的分辨率需要以下步骤：
import nctoolkit as nc
filename = '/Users/Nick/Desktop/SST/SST.nc'
data = nc.open_data(filename)
data.to_latlon(lon = [-179.5, 179.5], lat = [-89.5, 89.5], res = [1,1])
# visualize the data
data.plot()
关于你最初提到的 scipy.interpolate.griddata，我来给你解释一下：
仔细看看这个函数的参数说明（比如在SciPy文档里），确保你的输入数组有正确的形状。你可能需要做一些类似于：
import numpy as np
points = np.vstack([a.flat for a in np.meshgrid(lons,lats)]).T # (n,D)
values = sst.ravel() # (n)
等等。
我通常会用拉普拉斯滤波器来平滑我的数据。你可以试试下面这个函数，看看它是否对你的数据有帮助。这个函数可以选择带掩码或不带掩码来调用（比如针对海洋数据点的陆地/海洋掩码）。希望这对你有帮助。
# Laplace filter for 2D field with/without mask
# M = 1 on - cells used
# M = 0 off - grid cells not used
# Default is without masking
import numpy as np
def laplace_X(F,M):
jmax, imax = F.shape
# Add strips of land
F2 = np.zeros((jmax, imax+2), dtype=F.dtype)
F2[:, 1:-1] = F
M2 = np.zeros((jmax, imax+2), dtype=M.dtype)
M2[:, 1:-1] = M
MS = M2[:, 2:] + M2[:, :-2]
FS = F2[:, 2:]*M2[:, 2:] + F2[:, :-2]*M2[:, :-2]
return np.where(M > 0.5, (1-0.25*MS)*F + 0.25*FS, F)
def laplace_Y(F,M):
jmax, imax = F.shape
# Add strips of land
F2 = np.zeros((jmax+2, imax), dtype=F.dtype)
F2[1:-1, :] = F
M2 = np.zeros((jmax+2, imax), dtype=M.dtype)
M2[1:-1, :] = M
MS = M2[2:, :] + M2[:-2, :]
FS = F2[2:, :]*M2[2:, :] + F2[:-2, :]*M2[:-2, :]
return np.where(M > 0.5, (1-0.25*MS)*F + 0.25*FS, F)
# The mask may cause laplace_X and laplace_Y to not commute
# Take average of both directions
def laplace_filter(F, M=None):
if M == None:
M = np.ones_like(F)
return 0.5*(laplace_X(laplace_Y(F, M), M) +
laplace_Y(laplace_X(F, M), M))
如果你把数据用比如说双线性插值的方法重新调整到一个更粗的经纬度网格上，这样做会让数据看起来更平滑。
NCAR气候数据指南里有一个不错的重新调整网格的介绍（内容比较通用，不特指Python）。
在Python中，最强大的重新调整网格的工具是地球系统建模框架（ESMF）Python接口（ESMPy）。如果这个工具对你来说有点复杂，你可以看看
- EarthPy上关于重新调整网格的教程（比如使用Pyresample、cKDTree，或者Basemap）。
- 把你的数据转换成Iris立方体，并使用Iris的重新调整网格功能。
也许可以先看看使用Basemap的EarthPy重新调整网格教程，因为你已经在使用它了。
在你的例子中，做法是
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)