WRF如何转换投影+模拟台风路径可视化
先说结论
1、目前最推荐的方法是xesmf转换，插值方法多，自定义高，入门快 ，而且可导出权重文件进行复用
2、各种插值方式最常用的是线性或双线性插值，其他插值方法例如最邻近插值对边缘的处理一眼假，cubic慢
3、pyproj加scipy的griddata是第二推荐，进行pyproj投影转换后三种插值方法差别不明显，比之直接插值效果好
4、可视化仅作对比参考，现cartopy绘图能直接换投影
读取数据
import xarray as xr
ds = xr.open_dataset('/home/mw/input/wrf8852/sim/wrfout_d01_2019-08-08_19_00_00')
data = ds.T #取一层温度数据
data
1.1 仅使用griddata进行投影插值
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
XLONG=data.XLONG.data
XLAT =data.XLAT.data
datanum=data.data
# 输入wrf的网格
x = XLONG.flatten()
y = XLAT.flatten()
z = datanum.flatten()
# 定义目标网格
lon_target = np.arange(115, 130+1, 0.1)
lat_target = np.arange(22, 30+1, 0.1)
lon_target_grid, lat_target_grid = np.meshgrid(lon_target, lat_target)
# 进行网格插值
z_target_grid = griddata((x, y), z, (lon_target_grid, lat_target_grid), method='linear')
# 绘制投影后的数据
plt.pcolormesh(z_target_grid)
plt.colorbar()
plt.show()
1.2 多种griddata插值方法对比
In :
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
# 假设你已经有了以下数据
XLONG = data.XLONG.data
XLAT = data.XLAT.data
datanum = data.data
# 输入wrf的网格
x = XLONG.flatten()
y = XLAT.flatten()
z = datanum.flatten()
# 定义目标网格
lon_target = np.arange(115, 130+1, 0.1)
lat_target = np.arange(22, 30+1, 0.1)
lon_target_grid, lat_target_grid = np.meshgrid(lon_target, lat_target)
# 定义插值方法
methods = ['linear', 'nearest', 'cubic']
# 创建子图
fig, axs = plt.subplots(1, len(methods), figsize=(20, 5))
# 设置统一的色卡范围
vmin = min(z)
vmax = max(z)
# 进行网格插值并绘制子图
for i, method in enumerate(methods):
# 进行网格插值
z_target_grid = griddata((x, y), z, (lon_target_grid, lat_target_grid), method=method)
# 绘制子图，使用imshow来绘制颜色图
im = axs<i>.pcolormesh(z_target_grid, vmin=vmin, vmax=vmax)
axs<i>.set_title(method)
# 添加颜色条
fig.colorbar(im, ax=axs, orientation='vertical')
# 添加整体标题
plt.suptitle('Interpolation Methods')
# 显示图形
plt.show()</i></i>
'linear'：线性插值是一种基于线性关系进行插值的方法。它假设数据点之间的变化是线性的，并在相邻数据点之间进行插值。
'nearest'：最近邻插值是一种简单的插值方法，它将目标位置最近的数据点的值分配给目标位置。
'cubic'：三次插值是一种更复杂的插值方法，它基于数据点周围的局部曲线拟合进行插值。
这三种插值方法在速度、平滑度和准确性方面有所差异。通常情况下，'linear'插值速度较快，但在数据变化剧烈的地方可能会导致较大的误差；'nearest'插值计算速度快，但可能导致表面出现块状的不连续性；'cubic'插值在平滑度和准确性方面通常表现较好，但计算速度相对较慢。具体使用哪种插值方法应根据数据特点和需求进行选择
1.3 加入pyproj投影转换后griddata插值的多种方式对比
In :
import pyproj
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
# 输入wrf的网格
x = data.XLONG.data.flatten()
y = data.XLAT.data.flatten()
z = data.data.flatten()
# 定义Lambert投影和经纬度投影
wrf_proj = pyproj.Proj(proj='lcc', # projection type: Lambert Conformal Conic
lat_1=ds.TRUELAT1, lat_2=ds.TRUELAT2, # Cone intersects with the sphere
lat_0=ds.MOAD_CEN_LAT, lon_0=ds.STAND_LON, # Center point
a=6370000, b=6370000) # This is it! The Earth is a perfect sphere
wgs_proj = pyproj.Proj(proj='latlong', datum='WGS84')
# pyproj.transform() 函数用于将经纬度坐标(ds.CEN_LON, ds.CEN_LAT)从WGS84投影到WRF模型使用的投影坐标系。
e, n = pyproj.transform(wgs_proj, wrf_proj, ds.CEN_LON, ds.CEN_LAT)
dx, dy = ds.DX, ds.DY
nx, ny = ds.dims['west_east'], ds.dims['south_north']
# 通过计算网格的起始点（左下角）的坐标 x0 和 y0，基于网格的尺寸、分辨率和中心点坐标计算
x0 = -(nx-1) / 2. * dx + e
y0 = -(ny-1) / 2. * dy + n
# 用 np.meshgrid() 创建了一个二维网格 (xx, yy)，其中包含了整个模型的网格坐标信息
xx, yy = np.meshgrid(np.arange(nx) * dx + x0, np.arange(ny) * dy + y0)
# 使用 pyproj.transform() 将这些网格坐标点从 WRF 模型的投影坐标系转换回经纬度坐标系（PlateCarree投影），结果存储在 our_lons 和 our_lats
our_lons, our_lats = pyproj.transform(wrf_proj, wgs_proj, xx, yy)
# 定义插值方法
methods = ['linear', 'nearest', 'cubic']
# 创建子图
fig, axs = plt.subplots(1, len(methods), figsize=(20, 5))
# 进行网格插值并绘制子图
for i, method in enumerate(methods):
# 进行网格插值
z_target_grid = griddata((x, y), z, (our_lons, our_lats), method=method)
# 绘制子图
axs<i>.pcolormesh(our_lons, our_lats, z_target_grid)
axs<i>.set_title(method)
# 添加颜色条和整体标题
plt.suptitle('Interpolation Methods')
# 显示图形
plt.show()
</i></i>
2.1 xesmf转换投影
In :
import xesmf as xe
import xarray as xr
import meteva.base as meb
import numpy as np
XLAT = ds.XLAT.values
XLONG = ds.XLONG.values
# 输入wrf的网格
ds_in = xr.Dataset(
{
'lon': (['south_north', 'west_east'], XLONG),
'lat': (['south_north', 'west_east'], XLAT),
}
)
ds_out = xr.Dataset(
{
'lat': (['lat'], np.arange(22, 30+1, 0.1)),
'lon': (['lon'], np.arange(115, 130+1, 0.1)),
}
)
regridder = xe.Regridder(ds_in, ds_out, 'bilinear')
t_new = regridder(data)
# 绘制投影后的数据
plt.pcolormesh(t_new)
plt.colorbar()
plt.show()
2.2 使用xesmf的多种插值法进行转换对比
import xesmf as xe
import xarray as xr
import matplotlib.pyplot as plt
import numpy as np
XLAT = ds.XLAT.values
XLONG = ds.XLONG.values
# 输入wrf的网格
ds_in = xr.Dataset(
{
'lon': (['south_north', 'west_east'], XLONG),
'lat': (['south_north', 'west_east'], XLAT),
}
)
ds_out = xr.Dataset(
{
'lat': (['lat'], np.arange(22, 30+1, 0.1)),
'lon': (['lon'], np.arange(115, 130+1, 0.1)),
}
)
#conservative插值法使用时提示缺lon_b
methods = ['bilinear', 'patch','nearest_s2d','nearest_d2s']
# 创建子图
fig, axs = plt.subplots(1, len(methods), figsize=(30, 5))
# 进行插值并绘制子图
for i, method in enumerate(methods):
regridder = xe.Regridder(ds_in, ds_out, method)
t_new = regridder(data)
# 绘制子图
im = axs<i>.pcolormesh(t_new)
axs<i>.set_title(method)
# 添加颜色条
plt.colorbar(im, ax=axs<i>)
# 显示图形
plt.show()</i></i></i>
双线性插值（bilinear）：
优点：双线性插值是一种平滑的插值方法，在计算过程中考虑了周围四个网格点的权重。它在保持数据平滑性的同时，能够提供较为精确的插值结果。
缺点：尽管双线性插值是一种较为常用的插值方法，但在处理不规则或非均匀网格时可能会引入一些误差。
Patch插值（patch）：
优点：Patch 插值是一种多步骤的插值方法，通过将目标区域分成多个小块并进行插值，可以更好地处理不规则网格和不连续数据。它能够提供较高的插值精度。
缺点：Patch 插值需要较高的计算复杂度，并且可能需要较长的计算时间。
最近邻插值（nearest_s2d）：
优点：最近邻插值是一种简单快速的插值方法，它直接使用最近的一个源网格点的值来进行插值，不涉及其他计算。这种方法在处理离散数据或需要保留原始数据特征的情况下较为适用。
缺点：最近邻插值无法提供平滑的插值结果，可能导致插值值的不连续性，并且对于密集网格而言可能会引入一些误差。
反转最近邻插值（nearest_d2s）：
优点：反转最近邻插值是最近邻插值的一种变体，它根据目标网格上离坐标更近的源网格点来进行插值。这种方法可以在某种程度上避免最近邻插值带来的不连续性，并提供稍微平滑的插值结果。
缺点：反转最近邻插值在处理密集或高分辨率网格时可能会导致计算复杂度较高的问题，并且在插值过程中可能存在一定的误差。
扩展：利用WRF数据与台风实况表格进行模拟台风路径与实况对比
完整文件与代码在此
文章来源于微信公众号：气ython风雨
页:
[1]