关于WRF插值站点的二三事
前言
很多时候我们需要拿模拟数据和站点图作对比，那就需要把模拟数据插值到站点
今天来尝试两种WRF数据插值到站点的方法并使用meteva进行简单绘图
方法一：xesmf库重插值后使用meteva进行双线性插值到站点
方法二：proj+scipy重插值后使用meteva进行最临近插值到站点
import meteva.base as meb
import matplotlib.pyplot as plt
#由于meteva函数调用的是宋体，当前镜像的matplotlib字体库无宋体，先设置现有的tff
plt.rcParams['font.sans-serif'] = ['Source Han Sans CN']
plt.rcParams['axes.unicode_minus'] = False# 用来正常显示负号
station = meb.read_sta_alt_from_micaps3(meb.station_国家站)
方法一 xesmf库重插值
In :
import xarray as xr
wrffile = "/home/mw/input/wrfout3385/wrfout_d02_2022-07-14_0600.nc"
ds_wrf = xr.open_dataset(wrffile)
# 提取降雨量
temp = ds_wrf.Ttemp.plot()
Out:
<matplotlib.collections.QuadMesh at 0x7efda4ef71d0>
In :
import xesmf as xe
import xarray as xr
import meteva.base as meb
import numpy as np
XLAT = ds_wrf.XLAT.values
XLONG = ds_wrf.XLONG.values
# 输入wrf的网格
ds_in = xr.Dataset(
{
'lon': (['south_north', 'west_east'], XLONG),
'lat': (['south_north', 'west_east'], XLAT),
}
)
ds_in
In :
#设计输出网格
ds_out = xr.Dataset(
{
'lat': (['lat'], np.arange(27, 34+1, 0.1)),
'lon': (['lon'], np.arange(100, 111+1, 0.1)),
}
)
ds_out
In :#生成重插值网格
regridder = xe.Regridder(ds_in, ds_out, 'patch')
Create weight file: patch_90x90_80x120.nc
In :
#进行重插值
t_new = regridder(ds_wrf.T)
t_new.plot()
Out:
<matplotlib.collections.QuadMesh at 0x7efd9ae8e310>
数据重组装
In :
grid0 = meb.grid(,,gtime=["2022071406"],dtime_list = ,level_list = ,member_list = ["wrf"])
print(grid0)
tnewnew= meb.grid_data(grid0,t_new.data)
members:['wrf']
levels:
gtime:['20220714060000', '20220714060000', '1h']
dtimes:
glon:
glat:
meteva插值并可视化
In :
## 双线性插值
sta = meb.interp_gs_linear(tnewnew,station)
meb.tool.plot_tools.scatter_sta(sta)
方法二 pyproj+scipy重插值
In :
import pyproj
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
# 输入wrf的网格
x = temp.XLONG.data.flatten()
y = temp.XLAT.data.flatten()
z = temp.data.flatten()
# 定义Lambert投影和经纬度投影
wrf_proj = pyproj.Proj(proj='lcc', # projection type: Lambert Conformal Conic
lat_1=ds_wrf.TRUELAT1, lat_2=ds_wrf.TRUELAT2, # Cone intersects with the sphere
lat_0=ds_wrf.MOAD_CEN_LAT, lon_0=ds_wrf.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_wrf.CEN_LON, ds_wrf.CEN_LAT)
dx, dy = ds_wrf.DX, ds_wrf.DY
nx, ny = ds_wrf.dims['west_east'], ds_wrf.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)
# 进行网格插值
z_target_grid = griddata((x, y), z, (our_lons, our_lats), method='cubic')
# 绘制投影后的数据
plt.pcolormesh(z_target_grid)
plt.colorbar()
plt.show()
创建xarray数组
In :
# 创建xarray数据结构
t = xr.DataArray(z_target_grid,
coords=[('lat', our_lats[:, 0]), ('lon', our_lons)],
dims=['lat', 'lon'])
# 将DataArray添加到Dataset中
ds_inter = xr.Dataset({'temp': t})
ds_inter
meteva 转换xarray为grid_data（meteva可以绘制的格式）
In :
tnn =meb.xarray_to_griddata(ds_inter)
print(tnn) #对于da0里面的维度坐标名称为lon,lat，程序能够自动识别
插值可视化
In :
## 最临近插值
sta1 = meb.interp_gs_nearest(tnn,station)
meb.tool.plot_tools.scatter_sta(sta1)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-13-eab57681408a> in <module>
1 ## 最临近插值
2 sta1 = meb.interp_gs_nearest(tnn,station)
----> 3 meb.tool.plot_tools.scatter_sta(sta1)
/opt/conda/lib/python3.7/site-packages/meteva/base/tool/plot_tools.py in scatter_sta(sta0, value_column, map_extend, add_county_line, add_worldmap, clevs, cmap, fix_size, threshold, mean_value, print_max, print_min, save_dir, save_path, show, dpi, title, sup_fontsize, height, width, min_spot_value, grid, subplot, ncol, point_size, sup_title)
700 vmin_v = np.min(sta_without_iv.values)
701
--> 702 cmap1,clevs1 = meteva.base.tool.color_tools.def_cmap_clevs(cmap=cmap,clevs=clevs,vmin=vmin_v,vmax = vmax_v)
703 #clevs1, cmap1 = meteva.base.tool.color_tools.def_cmap_clevs(clevs=clevs, cmap=cmap, vmin = None, vmax=None)
704 #meteva.base.tool.color_tools.show_cmap_clev(cmap1,clevs1)
/opt/conda/lib/python3.7/site-packages/meteva/base/tool/color_tools.py in def_cmap_clevs(cmap, clevs, vmin, vmax, cut_colorbar)
958 vmax = vmin + 1.1
959 dif = (vmax - vmin) / 10.0
--> 960 inte = math.pow(10, math.floor(math.log10(dif)));
961 # 用基本间隔，将最大最小值除于间隔后小数点部分去除，最后把间隔也整数化
962 r = dif / inte
ValueError: cannot convert float NaN to integer出现nan值无法绘图，将nan值替换为0
In :sta1 = sta1.fillna(0)
In
meb.tool.plot_tools.scatter_sta(sta1)time or dtime or level 格式错误，请更改相应数据格式或直接指定title
以上可视化仅仅是展示插值后成果，需要进一步可视化可以使用matplotlib或者参考两种micaps站点数据的简单绘制方法
就使用而言，xesmf无疑是更简单的，并且插值后直接是xarray数组省去一步。
因为使用的插值方法不同就不作比较了，xesmf和griddata都有几种插值方法，感兴趣的读者可自行探索。
实际上在meteva的插值就使用了两种：最临近插值与双线性插值。效果好坏还需大家自行试验。
完整文件与代码在此
文章来源于微信公众号：气ython风雨
页:
[1]