重新网格化常规 netcdf 数据答案

【问题标题】：Regridding regular netcdf data重新网格化常规 netcdf 数据
【发布时间】：2023-03-26 21:03:01
【问题描述】：

我有一个包含全球海面温度的 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：条件必须是一维数组”。

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 的相对初学者，所以我不确定我做错了什么（或者可能有什么更好的方法）。任何建议表示赞赏！

【问题讨论】：

你有没有尝试使用conturf，我发现它可以产生更平滑的图
嗨，是的，我也尝试过 contourf，但它并不能真正解决问题。对于高分辨率数据，轮廓之间的边界并不平滑。如果我可以先对数据进行平滑处理，然后进行轮廓处理就可以了，但我还没有找到一个好的方法。
如何将精确值传递给contourf，在这些轮廓上必须绘制轮廓，间隔非常小，应该可以平滑地图上的数据。

标签： python numpy matplotlib matplotlib-basemap

【解决方案1】：

如果您在 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()

【讨论】：

【解决方案2】：

用 xarray 看这个例子... 使用ds.interp 方法并指定新的纬度和经度值。

http://xarray.pydata.org/en/stable/interpolation.html#example

【讨论】：

【解决方案3】：

也回答您关于scipy.interpolate.griddata 的原始问题：

仔细查看该函数的参数规范（例如，SciPy documentation）并确保您的输入数组具有正确的形状。您可能需要执行类似

的操作

import numpy as np
points = np.vstack([a.flat for a in np.meshgrid(lons,lats)]).T # (n,D)
values = sst.ravel() # (n)

等等

【讨论】：

【解决方案4】：

如果您将数据重新网格化 到更粗略的纬度/经度网格，例如使用双线性插值，这将产生一个平滑场。

NCAR ClimateData 指南有一个很好的introduction to regridding（通用，不是特定于 Python）。

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

EarthPy 重新网格化教程（例如，使用 Pyresample、cKDTree 或 Basemap）。
将您的数据转换为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)

【讨论】：

我发现scipy.interpolate.interp2d() 比底图快。
镇上的新人是xESMF - ESMF 接口到游戏规则改变者xarray。
xESMF 需要 Python3.5 或以上版本。

【解决方案5】：

我通常通过拉普拉斯过滤器运行我的数据以进行平滑处理。也许您可以尝试下面的功能，看看它是否对您的数据有帮助。可以使用或不使用掩码调用该函数（例如，海洋数据点的陆地/海洋掩码）。希望这可以帮助。 T

# 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))

【讨论】：