关于WRF插值站点的二三事

用户11172986

发布于 2024-06-20 16:45:36

760

发布于 2024-06-20 16:45:36

前言

很多时候我们需要拿模拟数据和站点图作对比，那就需要把模拟数据插值到站点今天来尝试两种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 [3]:

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.T

temp[0,0].plot()

Out[3]:

<matplotlib.collections.QuadMesh at 0x7efda4ef71d0>

In [4]:

import xesmf as xe
import xarray as xr
import meteva.base as meb
import numpy as np
XLAT = ds_wrf.XLAT.values[0, :, :]
XLONG = ds_wrf.XLONG.values[0, :, :]
# 输入wrf的网格
ds_in = xr.Dataset(  
                    {
                     'lon': (['south_north', 'west_east'], XLONG),
                     'lat': (['south_north', 'west_east'], XLAT), 
                    }
                          )
ds_in

In [5]:

#设计输出网格
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 [17]:

#生成重插值网格
regridder = xe.Regridder(ds_in, ds_out, 'patch')

Create weight file: patch_90x90_80x120.nc

In [18]:

#进行重插值
t_new = regridder(ds_wrf.T[0,0])
t_new.plot()

Out[18]:

<matplotlib.collections.QuadMesh at 0x7efd9ae8e310>

数据重组装

In [19]:

grid0 = meb.grid([100,111.9,0.1],[27,34.9,0.1],gtime=["2022071406"],dtime_list = [0],level_list = [0],member_list = ["wrf"])
print(grid0)
tnewnew= meb.grid_data(grid0,t_new.data)

members:['wrf']
levels:[0]
gtime:['20220714060000', '20220714060000', '1h']
dtimes:[0]
glon:[100, 111.9, 0.1]
glat:[27, 34.9, 0.1]

meteva插值并可视化

In [20]:

## 双线性插值
sta = meb.interp_gs_linear(tnewnew,station)
meb.tool.plot_tools.scatter_sta(sta)

方法二 pyproj+scipy重插值

In [10]:

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[0,3].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 [11]:

# 创建xarray数据结构
t = xr.DataArray(z_target_grid,
                  coords=[('lat', our_lats[:, 0]), ('lon', our_lons[0, :])],
                  dims=['lat', 'lon'])

# 将DataArray添加到Dataset中
ds_inter = xr.Dataset({'temp': t})
ds_inter

meteva 转换xarray为grid_data（meteva可以绘制的格式）

In [12]:

tnn =meb.xarray_to_griddata(ds_inter)
print(tnn) #对于da0里面的维度坐标名称为lon,lat，程序能够自动识别

***
<xarray.DataArray 'data0' (member: 1, level: 1, time: 1, dtime: 1, lat: 91, lon: 90)>
array([[[[[[        nan,         nan, 20.69002073, ...,  8.90399785,
             8.79607043,  8.73185356],
           [        nan,         nan, 21.0672533 , ...,  8.95046101,
             8.8558546 ,  8.7926331 ],
           [        nan,         nan, 21.44727255, ...,  8.99169476,
             8.89961342,  8.83765678],
           ...,
           [        nan,         nan, 25.85778439, ..., 11.51698046,
            11.63704809, 11.71500662],
           [        nan,         nan,         nan, ...,         nan,
                    nan,         nan],
           [        nan,         nan,         nan, ...,         nan,
                    nan,         nan]]]]]])
Coordinates:
  * member   (member) int64 0
  * level    (level) int64 0
  * time     (time) datetime64[ns] 1970-01-01
  * dtime    (dtime) int64 0
  * lat      (lat) float64 27.24 27.33 27.42 27.51 ... 34.96 35.05 35.13 35.22
  * lon      (lon) float64 100.9 101.0 101.1 101.2 ... 109.5 109.6 109.7 109.8

插值可视化

In [13]:

## 最临近插值
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[plot_data_names].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 [14]:

sta1 = sta1.fillna(0)

In [15]:

meb.tool.plot_tools.scatter_sta(sta1)

time or dtime or level 格式错误，请更改相应数据格式或直接指定title

以上可视化仅仅是展示插值后成果，需要进一步可视化可以使用matplotlib或者参考两种micaps站点数据的简单绘制方法就使用而言，xesmf无疑是更简单的，并且插值后直接是xarray数组省去一步。因为使用的插值方法不同就不作比较了，xesmf和griddata都有几种插值方法，感兴趣的读者可自行探索。实际上在meteva的插值就使用了两种：最临近插值与双线性插值。效果好坏还需大家自行试验。

完整文件与代码在此

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2023-12-09，如有侵权请联系 cloudcommunity@tencent.com 删除

可视化