前言
博主在早期对meteva的使用写了一个笔记，就是meteva，这可能是气象萌新最需要的python库
在使用中发现它不能对有兰伯特投影的wrfout数据直接绘图，所以使用了其他库进行重新网格插值再绘图
今天在逛meteva的showdoc时刷新出了一个官方教程，大体是将wrfout数据转为pandas格式
然后使用idw进行插值绘图
下面让我们开始实践吧
温馨提示
由于可视化代码过长隐藏，可点击回旋镖！meteva也能绘制wrfout气象要素分布运行Fork查看
🔜🔜若没有成功加载可视化图，点击运行可以查看
ps:隐藏代码在【代码已被隐藏】所在行，点击所在行，可以看到该行的最右角，会出现个三角形，点击查看即可
导入库与读取数据
In [23]:
import xarray as xrimport matplotlib.pyplot as plt#由于meteva函数调用的是宋体，当前镜像的matplotlib字体库无宋体，先设置现有的tffplt.rcParams['font.sans-serif'] = ['Source Han Sans CN']plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号path = "/home/mw/input/wrfout3385/wrfout_d02_2022-07-14_1100.nc"grd0 = xr.open_dataset(path) #通过xarray读入网格数据print(grd0)
<xarray.Dataset> Size: 51MB
Dimensions: (Time: 1, south_north: 90, west_east: 90,
bottom_top: 49, bottom_top_stag: 50,
soil_layers_stag: 4, west_east_stag: 91,
south_north_stag: 91, seed_dim_stag: 8)
Coordinates:
XLAT (Time, south_north, west_east) float32 32kB ...
XLONG (Time, south_north, west_east) float32 32kB ...
XTIME (Time) datetime64[ns] 8B ...
XLAT_U (Time, south_north, west_east_stag) float32 33kB ...
XLONG_U (Time, south_north, west_east_stag) float32 33kB ...
XLAT_V (Time, south_north_stag, west_east) float32 33kB ...
XLONG_V (Time, south_north_stag, west_east) float32 33kB ...
Dimensions without coordinates: Time, south_north, west_east, bottom_top,
bottom_top_stag, soil_layers_stag,
west_east_stag, south_north_stag, seed_dim_stag
Data variables: (12/213)
Times (Time) |S19 19B ...
LU_INDEX (Time, south_north, west_east) float32 32kB ...
ZNU (Time, bottom_top) float32 196B ...
ZNW (Time, bottom_top_stag) float32 200B ...
ZS (Time, soil_layers_stag) float32 16B ...
DZS (Time, soil_layers_stag) float32 16B ...
... ...
PCB (Time, south_north, west_east) float32 32kB ...
PC (Time, south_north, west_east) float32 32kB ...
LANDMASK (Time, south_north, west_east) float32 32kB ...
LAKEMASK (Time, south_north, west_east) float32 32kB ...
SST (Time, south_north, west_east) float32 32kB ...
SST_INPUT (Time, south_north, west_east) float32 32kB ...
Attributes: (12/134)
TITLE: OUTPUT FROM WRF V4.4.2 MODEL
START_DATE: 2022-07-14_00:00:00
SIMULATION_START_DATE: 2022-07-14_00:00:00
WEST-EAST_GRID_DIMENSION: 91
SOUTH-NORTH_GRID_DIMENSION: 91
BOTTOM-TOP_GRID_DIMENSION: 50
... ...
ISLAKE: 21
ISICE: 15
ISURBAN: 13
ISOILWATER: 14
HYBRID_OPT: 2
ETAC: 0.2
转为pandas
In [25]:
import pandas as pdimport numpy as npsta = pd.DataFrame({"level":0, "time":grd0.coords["XTIME"].values[0],"dtime":0,"id":np.arange(grd0.coords["XLAT"].values.size),"lon":grd0.coords["XLONG"].values.flatten(),"lat":grd0.coords["XLAT"].values.flatten(),"RAIN":grd0.variables["RAINC"].values.flatten()+grd0.variables["RAINNC"].values.flatten()})
print(sta)
level time dtime id lon lat RAIN
0 0 2022-07-14 11:00:00 0 0 100.901947 27.242496 0.323906
1 0 2022-07-14 11:00:00 0 1 101.002350 27.238125 0.238994
2 0 2022-07-14 11:00:00 0 2 101.102722 27.233673 0.315282
3 0 2022-07-14 11:00:00 0 3 101.203064 27.229111 0.267089
4 0 2022-07-14 11:00:00 0 4 101.303406 27.224480 0.215387
... ... ... ... ... ... ... ...
8095 0 2022-07-14 11:00:00 0 8095 110.649567 34.472286 19.087238
8096 0 2022-07-14 11:00:00 0 8096 110.757538 34.459393 17.735456
8097 0 2022-07-14 11:00:00 0 8097 110.865448 34.446400 15.124860
8098 0 2022-07-14 11:00:00 0 8098 110.973358 34.433327 13.774797
8099 0 2022-07-14 11:00:00 0 8099 111.081207 34.420151 13.101636
[8100 rows x 7 columns]
IDW插值
In [26]:
import meteva.base as mebgrid = meb.grid([100,112,0.05],[28,36,0.05]) # 设置一个等经纬度的网格信息变量grd1 = meb.interp_sg_idw(sta,grid = grid,nearNum=4,effectR=20) # 将站点数据插值到网格上print(grd1)
<xarray.DataArray 'data0' (member: 1, level: 1, time: 1, dtime: 1, lat: 161,
lon: 241)> Size: 310kB
array([[[[[[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]]]]]])
Coordinates:
* member (member) <U4 16B 'RAIN'
* level (level) int64 8B 0
* time (time) datetime64[ns] 8B 2022-07-14T11:00:00
* dtime (dtime) int64 8B 0
* lat (lat) float64 1kB 28.0 28.05 28.1 28.15 ... 35.85 35.9 35.95 36.0
* lon (lon) float64 2kB 100.0 100.0 100.1 100.2 ... 111.9 112.0 112.0
Attributes:
data_start_columns: 6
绘图
In [27]:
import cmapscmaps0=cmaps.prcp_1meb.contourf_2d_grid(grd1,cmap=cmaps0) #绘制插值结果
高度层绘制
按照官方的思路，我们对高度层的混合比数据进行绘制
In [28]:
from netCDF4 import Datasetimport numpy as npimport pandas as pdfrom wrf import (to_np, getvar, ALL_TIMES, CoordPair, interplevel)#由于meteva函数调用的是宋体，当前镜像的matplotlib字体库无宋体，先设置现有的tffplt.rcParams['font.sans-serif'] = ['Source Han Sans CN']plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号wrf_file = '/home/mw/input/wrfout3385/wrfout_d02_2022-07-14_0600.nc'ncfile = Dataset(wrf_file)# 使用wrf-python的getvar函数获取整个时间序列的湿度数据qvapor = getvar(ncfile, 'QVAPOR', timeidx=ALL_TIMES, method='cat')# 定义需要插值的目标气压值，这里为500hPatarget_plev = 500.0 # in hPa# 获取模型中的气压数据，这是进行插值所必需的pressure = getvar(ncfile, 'P')# 使用interplevel函数进行垂直插值，得到500hpa高度的QVAPOR数据qvapor_500hpa = interplevel(qvapor, pressure, target_plev)# 获取必要的坐标信息time_coord = ncfile.START_DATElat = getvar(ncfile, 'XLAT').valueslon = getvar(ncfile, 'XLONG').values# 将数据整理到Pandas DataFramesta1 = pd.DataFrame({"level": target_plev, # 直接指定为500hPa"time": [time_coord][0], "dtime":0,"id": np.arange(lat.size).flatten(), # 网格点ID，可选"lon": lon.flatten(),"lat": lat.flatten(),"QVAPOR": qvapor_500hpa.values.flatten(), # 使用插值后的QVAPOR作为湿度数据})print(sta1.head()) # 打印DataFrame的前几行以检查结果# 关闭NetCDF文件
ncfile.close()
level time dtime id lon lat QVAPOR
0 500.0 2022-07-14_00:00:00 0 0 100.901947 27.242496 0.001734
1 500.0 2022-07-14_00:00:00 0 1 101.002350 27.238125 0.001666
2 500.0 2022-07-14_00:00:00 0 2 101.102722 27.233673 0.001562
3 500.0 2022-07-14_00:00:00 0 3 101.203064 27.229111 0.001444
4 500.0 2022-07-14_00:00:00 0 4 101.303406 27.224480 0.001443
In [29]:
grid = meb.grid([100,112,0.05],[28,35,0.05]) # 设置一个等经纬度的网格信息变量grd1 = meb.interp_sg_idw(sta1,grid = grid,nearNum=4,effectR=20) # 将站点数据插值到网格上
print(grd1)
<xarray.DataArray 'data0' (member: 1, level: 1, time: 1, dtime: 1, lat: 141,
lon: 241)> Size: 272kB
array([[[[[[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]]]]]])
Coordinates:
* member (member) <U6 24B 'QVAPOR'
* level (level) int64 8B 500
* time (time) datetime64[ns] 8B 2022-07-14
* dtime (dtime) int64 8B 0
* lat (lat) float64 1kB 28.0 28.05 28.1 28.15 ... 34.85 34.9 34.95 35.0
* lon (lon) float64 2kB 100.0 100.0 100.1 100.2 ... 111.9 112.0 112.0
Attributes:
data_start_columns: 6
下面我绘制了基于cartopy的两种版本的500hPa混合比分布，可作对比
等经纬度网格版本
兰伯特版本
小结
能摸索出新的方法对数据处理绘图是愉快的
图片空白部分不是数据有误或者作图有错，而是地形的原因
总的来说meteva绘图一句话，但是前期的数据转换较为复杂 ，还是会劝退萌新
