Python气象数据处理与绘图：4D 插值

1、导入模块

1import cartopy.crs
2import matplotlib
3import matplotlib.pyplot
4import numpy
5import pyinterp
6import pyinterp.backends.xarray
7import xarray
8import warnings
9warnings.filterwarnings('ignore')

2、读取数据

1ds = xarray.open_dataset("pres_temp_4D.nc")
2interpolator = pyinterp.backends.xarray.Grid4D(ds.pressure)

3、创建网格

1mx, my, mz, mu = numpy.meshgrid(numpy.arange(-125, -70, 0.5),
2                                numpy.arange(25, 50, 0.5),
3                                numpy.datetime64("2000-01-01T12:00"),
4                                0.5,
5                                indexing="ij")

4、选择插值方法（quadrivariate）

1quadrivariate = interpolator.quadrivariate(
2    dict(longitude=mx.flatten(),
3         latitude=my.flatten(),
4         time=mz.flatten(),
5         level=mu.flatten())).reshape(mx.shape)

5、选择插值方法（bicubic）

1interpolator = pyinterp.backends.xarray.Grid4D(ds.pressure, increasing_axes=True)
2bicubic = interpolator.bicubic(dict(longitude=mx.flatten(),
3                                    latitude=my.flatten(),
4                                    time=mz.flatten(),
5                                    level=mu.flatten()),
6                               nx=2,
7                               ny=2).reshape(mx.shape)

6、进行数据插值

1quadrivariate = quadrivariate.squeeze(axis=(2, 3))
2bicubic = bicubic.squeeze(axis=(2, 3))
3lons = mx[:, 0].squeeze()
4lats = my[0, :].squeeze()

7、可视化插值结果

 1fig = matplotlib.pyplot.figure(figsize=(5, 4))
 2ax1 = fig.add_subplot(211, projection=cartopy.crs.PlateCarree(central_longitude=180))
 3pcm = ax1.pcolormesh(lons,
 4                     lats,
 5                     quadrivariate.T,
 6                     cmap='jet',
 7                     transform=cartopy.crs.PlateCarree())
 8ax1.coastlines()
 9ax1.set_title("Trilinear")
10
11ax2 = fig.add_subplot(212, projection=cartopy.crs.PlateCarree(central_longitude=180))
12pcm = ax2.pcolormesh(lons,
13                     lats,
14                     bicubic.T,
15                     cmap='jet',
16                     transform=cartopy.crs.PlateCarree())
17ax2.coastlines()
18ax2.set_title("Spline & Linear in time")
19fig.colorbar(pcm, ax=[ax1, ax2], shrink=0.8)
20fig.show()