meteva，这可能是气象萌新最需要的python库

小葵花妈妈课堂开课了

本文内容参考meteva官方文档：                                                           https://www.showdoc.com.cn/meteva/
面向群体：拿着micaps数据却不知如何读取与绘图的小白
初衷：向广大气象学子介绍meteva库的用途，达到降低micaps门槛的目的，主打“我淋过雨我撑伞”
应用场景：比如说组会前赶图（

前言

许多同学在面对micaps、nc、grib等格式的气象数据傻眼，不知用什么方式读取他们，可视化他们。

老师也不可能一步步喂食，自身摸索耗时过多且盲目。

那么有没有一个既代码两少又能满足气象数据读取与可视化需求的库呢？

meteva

• 只需简单几行代码，你就能够读取、处理和操作各种格式的气象数据。

• 不管你是要进行数据筛选、插值、格点场计算还是统计分析，meteva库都提供了一系列酷炫的方法和工具，满足你的各种需求。

• meteva库有着强大的可视化功能，可以制作各种炫酷的气象图，比如等值线图、填色图、风场图等等。比起matplotlib动不动数十行的绘图设置，meteva表示一句搞定。

总而言之，meteva库是你的气象数据好伙伴！

它不仅能够让你轻松处理和可视化各种格式的数据，还能帮助你展现出充满活力的数据分析结果。

相信我，选择meteva库绝对能让你在气象数据处理的世界里燃起激情的火焰！Let's go，一起点燃数据的热情吧！

当新手开始使用meteva库时，可能会遇到一些常见的问题。下面我分为几个方面：

数据读取

如何读取micaps，nc，grib格式的气象数据？

案例：1.1 站点数据读取：以micaps3类数据为例
1.3 格点数据读取：以era5的nc数据为例
1.7 grib 数据处理简单流程：使用era5 girb格式为例

数据操纵

如何进行气象数据的插值操作？

案例 ：2.1 测试wrfout可视化：使用xesmf重插值后使用meteva可视化wrfout数据

案例 ：2.2 nc网格数据插值站点三步走：meteva读取nc数据并插值到气象站点

数据可视化

1. 如何绘制站点数据图？案例：1.2 站点数据绘制：以micaps3类数据为例

2. 如何绘制格点数据图？案例：1.4 格点数据绘制 ：以era5的nc数据为例

3. 地图白化如何操作？案例：1.5 地图白化

4. 多图层的气象要素要怎么实现？案例：1.6 多图层绘制：使用micaps3类与micaps4类数据为例

这些问题是新手在使用meteva库时可能遇到的一些疑惑。

来吧，一起在阅读和运行中寻找答案。

希望本文能给你微小的帮助。

METEVA

是什么

Meteva是一个纯python程序库，提供了常用的各种气象预报检验评估算法函数，气象检验分析的图片和表格型产品的制作函数，以及可以实用的检验评估系统代码示例。共包含检验基础算法100+项，检验应用工具20+项，并支持30+种分类方式开展精细化检验。

可以干什么

博主主要用于读取micaps数据与绘图 ,就业之后可用于业务评分。

扩展阅读

micaps的产品种类可以查看网站                                                           https://lyqx.sinaapp.com/wfiles/micaps-dir/

1.1 站点数据读取：以micaps3类数据为例读取地表温度数据

In [1]:

import meteva.base as meb

#读取地表2m温度filename = '/home/mw/input/meteva2260/t.000'#注意下面是micaps 3类数据station = meb.read_stadata_from_micaps3(filename)station.head()

可以看出meteva将micaps文档转为df表格，可大大方便了我们使用pandas处理

station.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 2289 entries, 0 to 2288
Data columns (total 7 columns):
 #   Column  Non-Null Count  Dtype         
---  ------  --------------  -----         
 0   level   2289 non-null   int64         
 1   time    2289 non-null   datetime64[ns]
 2   dtime   2289 non-null   int64         
 3   id      2289 non-null   int32         
 4   lon     2289 non-null   float64       
 5   lat     2289 non-null   float64       
 6   data0   2289 non-null   float64       
dtypes: datetime64[ns](1), float64(3), int32(1), int64(2)
memory usage: 116.4 KB

修改变量列名

In [5]:

#该函数是专门针对变量列的meb.set_stadata_names(station,["T2"])print(station)

      level                time  dtime     id     lon    lat    T2
0         0 2022-09-04 17:00:00      0  56473  102.77  28.95  16.8
1         0 2022-09-04 17:00:00      0  56263  101.88  30.88  11.9
2         0 2022-09-04 17:00:00      0  56593  104.92  28.58  23.1
3         0 2022-09-04 17:00:00      0  56491  104.57  28.70  21.7
4         0 2022-09-04 17:00:00      0  56290  104.18  30.78  19.6
...     ...                 ...    ...    ...     ...    ...   ...
2284      0 2022-09-04 17:00:00      0  58205  115.43  32.47  23.5
2285      0 2022-09-04 17:00:00      0  58208  115.62  32.17  22.0
2286      0 2022-09-04 17:00:00      0  58104  115.07  33.40  20.4
2287      0 2022-09-04 17:00:00      0  57294  114.02  32.83  20.2
2288      0 2022-09-04 17:00:00      0  53492  114.15  40.10  15.9

[2289 rows x 7 columns]

当读取产品在micaps的列号时需要查表可用下面

print(meb.m1_element_column.温度)print(meb.m1_element_column.气压)

19
8

1.2 站点数据绘图：以micaps3类数据为例绘制地表温度站点分布

#导入必要包import numpy as npimport matplotlib.pyplot as pltfrom meteva.base.tool.plot_tools import add_china_map_2basemap#由于meteva函数调用的是宋体，当前镜像的matplotlib字体库无宋体，先设置现有的tffplt.rcParams['font.sans-serif'] = ['Source Han Sans CN']  plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号

# 使用 scatter_sta 绘制站点分布meb.tool.plot_tools.scatter_sta(station)

meb.tool.plot_tools.scatter_sta(station,value_column=0,map_extend = [110,130,20,40], #设置范围                                cmap =meb.cmaps.temp_2m,  #设置色卡                             fix_size=True) #根据打点值的大小决定半径

meb.tool.plot_tools.scatter_sta(station,value_column=0,map_extend = [110,130,20,40], #设置范围                                cmap =meb.cmaps.temp_2m,  #设置色卡                             fix_size=True, #根据打点值的大小决定半径                             add_county_line=True) # 县界

1.3 格点数据读取：以era5的nc数据为例读取地表2m温度分布数据

grid0 = meb.grid([110,130,1],[20,40,1])  #设置范围应小于或等于era5数据的原有范围grd = meb.io.read_griddata_from_nc('/home/mw/input/meteva2260/air2m20210905.nc',grid = grid0)  #读取网格数据，用于测试meb.set_griddata_coords(grd,name = "Temper_2m",dtime_list = [24],member_list = ["ECMWF"])    #设置数据的时空属性和要素名称grd = grd -273.15  #因单位是开尔文，转为摄氏print(grd)

<xarray.DataArray 'Temper_2m' (member: 1, level: 1, time: 24, dtime: 1, lat: 21, lon: 21)>
array([[[[[[26.99346313, 26.90334473, 26.82338867, ..., 27.49712524,
            27.02071533, 27.84508667],
           [28.07521973, 27.12423096, 26.76287231, ..., 27.34880981,
            27.77855835, 28.21383057],
           [26.78829346, 26.84465942, 27.01101074, ..., 27.53411255,
            28.52852783, 28.72167358],
           ...,
           [16.88796387, 16.69436035, 14.26885376, ..., 19.58706055,
            21.86425171, 23.13405762],
           [16.23919067, 17.37825928, 14.09462891, ..., 20.38601074,
            22.11422119, 22.62807617],
           [14.19075928, 15.96010742, 13.08495483, ..., 20.96224365,
            21.48162231, 22.00518188]]],


         [[[27.47219238, 27.33168945, 26.93706665, ..., 27.86312256,
            27.32753906, 28.11495361],
           [28.50725708, 26.91951904, 26.84142456, ..., 27.87096558,
            27.51516113, 28.28869019],
           [27.25960693, 27.02303467, 27.47817383, ..., 27.94906006,
...
            22.10036621, 22.5102478 ],
           [13.55120239, 14.24480591, 12.09240112, ..., 20.92849121,
            22.21172485, 22.30321655],
           [11.11968384, 13.25640259, 12.18804321, ..., 19.87957153,
            21.51442871, 21.90490112]]],


         [[[27.24807129, 28.05440674, 28.85335693, ..., 28.62508545,
            28.56826172, 28.53682861],
           [27.13485107, 28.41991577, 28.42546997, ..., 28.67962036,
            28.59182129, 28.40792236],
           [26.34191284, 26.11687622, 26.21297607, ..., 28.58950195,
            28.61724243, 28.46105347],
           ...,
           [15.09551392, 14.96105347, 14.15426025, ..., 19.19890747,
            22.23623047, 22.69136963],
           [13.08773193, 14.08168945, 11.93944702, ..., 21.08791504,
            22.24593506, 22.42382202],
           [10.87108765, 12.66216431, 11.99996338, ..., 20.18361816,
            21.60732422, 21.91182861]]]]]])
Coordinates:
  * member   (member) <U5 'ECMWF'
  * level    (level) int64 0
  * time     (time) datetime64[ns] 2021-09-05 ... 2021-09-05T23:00:00
  * dtime    (dtime) int64 24
  * lat      (lat) int64 20 21 22 23 24 25 26 27 28 ... 33 34 35 36 37 38 39 40
  * lon      (lon) int64 110 111 112 113 114 115 116 ... 125 126 127 128 129 130

1.4 格点数据绘制：以era5的nc数据为例绘制地表2m温度分布数据

subplot: 当sta0包含了多个level,time,dtime或member的数据时，例如当subplot = level 时，则一张图中会包含多个层次数据的子图，而其它维度如果size!=1 则会以不同文件的形式依次生成

meb.tool.plot_tools.contourf_2d_grid(grd,                                     cmap =meb.cmaps.temp_2m, #色卡设置                                     subplot = "time", #按时间排序                                     ncol = 3, #设置三列                                     sup_fontsize = 9,#标题字体                                     sup_title = "测试") #标题

findfont: Font family ['Times New Roman'] not found. Falling back to DejaVu Sans.

1.5 地图白化：在1.4可视化基础上的白化操作

# 白化后显示广东与浙江的T2分布meb.tool.plot_tools.contourf_2d_grid(    grd, cmap=meb.cmaps.temp_2m, clip=[   #白化函数，选取区域        "guangdong", "zhejiang"], subplot="time", ncol=3,        sup_fontsize = 9)#标题字体

1.6 多图层绘制：使用micaps3类与micaps4类数据为例可视化相对湿度、海平面气压等值线、24小时降水站点分布

1.6.1 底图绘制 creat_axs

nplot: 子图的数量。
map_extend: 地图范围，格式为 [(left, right, bottom, top)]。
ncol: 列数，默认为 None，表示根据子图数量自动确定列数。
height: 子图的高度，默认为 None，表示自动计算高度。
width: 子图的宽度，默认为 None，表示自动计算宽度。
dpi: 图像的分辨率，默认为 300。
sup_title: 整个图像的标题，默认为 None。
sup_fontsize: 整个图像的标题字体大小，默认为 12。
add_county_line: 是否添加县界线，默认为 False。
add_worldmap: 是否添加世界地图，默认为 True。
title_list: 子图的标题列表，默认为 None。
add_index: 是否在子图上添加索引，默认为 None。

map_extend = [100,130,20,40]axs = meb.creat_axs(4,map_extend,ncol = 2,sup_title = "测试底图",add_index = ["a","b","c","d"],                    sup_fontsize = 8)

1.6.2 叠加站点数值 add_scatter_text

map_extend = [110, 120, 20, 25]axs = meb.creat_axs(    1,    map_extend,    sup_title="2022年9月4日t2",    add_index=[        "a"],    sup_fontsize=8)image = meb.add_scatter_text(axs[0], station, tag=0, font_size=5,cmap =meb.cmaps.temp_2m)#tag为数值转文字保留的有效位数

#读取海平面气压slp = meb.read_griddata_from_micaps4("/home/mw/input/meteva2260/slp.000")slp = meb.smooth(slp,1) #平滑

1.6.3 叠加散点图函数 meb.add_scatter

ax: Matplotlib的Axes对象，表示要在其上添加散点图。
map_extend: 地图范围，格式为[(left, right, bottom, top)]，用于设置坐标轴范围。
sta0: 站点数据，包含经纬度、值等信息。
cmap: 散点图的颜色映射，默认为None，表示使用默认的颜色映射。
clevs: 颜色刻度值，默认为None，表示自动计算刻度值。
point_size: 散点的大小，默认为None，表示使用默认大小。
fix_size: 是否使用固定大小的散点，默认为True。
title: 图像的标题，默认为None，表示无标题。
threshold: 阈值，用于显示大于阈值的散点，默认为2。
min_spot_value: 最小的散点值，默认为0。
mean_value: 平均值，用于计算高于平均值的散点，默认为2。
grid: 是否显示网格线，默认为False。
add_colorbar: 是否添加颜色条，默认为True。
alpha: 散点的透明度，默认为None，表示使用默认透明度。

rh500 = meb.read_griddata_from_micaps4('/home/mw/input/meteva2260/rh.000')map_extend = [100,130,20,40]axs = meb.creat_axs(2,map_extend,ncol = 2,sup_title = "2022年9月4日",                    add_index = ["a","b"],sup_fontsize = 8)image = meb.add_scatter(axs[0],map_extend,station,cmap = meb.cmaps.temp_2m,add_colorbar=True,alpha=1)image = meb.add_contour(axs[0],slp,clevs = np.arange(950,1020,4))image = meb.add_contour(axs[1],slp,clevs = np.arange(950,1020,4))

[ 950  954  958  962  966  970  974  978  982  986  990  994  998 1002
 1006 1010 1014 1018]
[ 950  954  958  962  966  970  974  978  982  986  990  994  998 1002
 1006 1010 1014 1018]

1.6.4 叠加马赛克图层 add_mesh

ax: matplotlib 的 Axes 对象，表示需要绘制等值线/填色图的子图。
grd: meteva 的 GridData 对象，表示网格数据。
cmap: 颜色映射表，默认为 'rainbow'。
clevs: 等值线的颜色分级，默认为 None，表示根据数据自动确定分级。
add_colorbar: 是否添加颜色条，默认为 True。

#读取24小时降水数据rain24 = meb.read_stadata_from_micaps3('/home/mw/input/meteva2260/ob_rain24_BT19123108.000')map_extend = [70,140,15,55]axs = meb.creat_axs(2,map_extend,ncol = 2,sup_title = "rain and rh500 test",                    add_index = ["a","b"],width=9,sup_fontsize = 12)image = meb.add_scatter(axs[0],map_extend,rain24,cmap = meb.cmaps.rain_24h,add_colorbar=True,alpha=1)image = meb.add_mesh(axs[1],rh500,cmap = meb.cmaps.rh,clevs = np.arange(0,101,10),add_colorbar=False)

findfont: Font family ['Times New Roman'] not found. Falling back to DejaVu Sans.

练习：自行学会使用叠加填色图层函数add_contourf与叠加风羽函数add_barbs。相关用法例子官网可查。

1.7 grib 数据处理简单流程：使用era5 girb格式为例可视化降水数据

In [19]:

#打印一个era5的grib数据看看数据结构import meteva.base as mebmeb.print_grib_file_info('/home/mw/input/meteva2260/ERA5202007130000_prep.grib')

/home/mw/input/meteva2260/ERA5202007130000_prep.grib中只有一种leve_type，
请根据以下数据内容信息，确认其中的level维度名称
<xarray.Dataset>
Dimensions:     (latitude: 721, longitude: 1440)
Coordinates:
    number      int64 ...
    time        datetime64[ns] ...
    step        timedelta64[ns] ...
    surface     float64 ...
  * latitude    (latitude) float64 90.0 89.75 89.5 89.25 ... -89.5 -89.75 -90.0
  * longitude   (longitude) float64 0.0 0.25 0.5 0.75 ... 359.2 359.5 359.8
    valid_time  datetime64[ns] ...
Data variables:
    tp          (latitude, longitude) float32 ...
Attributes:
    GRIB_edition:            1
    GRIB_centre:             ecmf
    GRIB_centreDescription:  European Centre for Medium-Range Weather Forecasts
    GRIB_subCentre:          0
    Conventions:             CF-1.7
    institution:             European Centre for Medium-Range Weather Forecasts
    history:                 2023-07-17T04:42 GRIB to CDM+CF via cfgrib-0.9.9...

#试图直接导入路径进行读取pre=meb.read_griddata_from_grib('/home/mw/project/ERA5202007130000_prep.grib')pre

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-20-a25a306c3ecb> in <module>
      1 #试图直接导入路径进行读取
----> 2 pre=meb.read_griddata_from_grib('/home/mw/project/ERA5202007130000_prep.grib')
      3 pre

TypeError: read_griddata_from_grib() missing 1 required positional argument: 'level_type'

提示缺少level_type，使用help()

#官网文档查到有些数据需要传入level

help(meb.read_griddata_from_grib)

Help on function read_griddata_from_grib in module meteva.base.io.read_griddata:

read_griddata_from_grib(filename, level_type, grid=None, value_name=None, member_dim=None, time_dim=None, dtime_dim=None, lat_dim=None, lon_dim=None, level=None, time=None, dtime=None, data_name='data0', filter_by_keys={}, show=False)

#level_type估计是填入地面或高空，填入surface试试pre=meb.read_griddata_from_grib('/home/mw/project/ERA5202007130000_prep.grib',level_type='surface')

#简单可视化meb.tool.plot_tools.contourf_2d_grid(pre)

1.8 绘图后如何保存图片以及指定格式保存

指定文件夹保存

meb.tool.plot_tools.scatter_sta(station,save_dir='/home/mw/project')

图片已保存至/home/mw/project/data0.png

指定文件夹并存为pdf格式

meb.tool.plot_tools.scatter_sta(station,save_path='/home/mw/project/data0.pdf')

图片已保存至/home/mw/project/data0.pdf

2.1 测试wrfout可视化：使用xesmf重插值后使用meteva可视化wrfout温度数据

直接使用meteva读取是不行的，当前不支持兰伯特投影

2.1.1 读取wrfout数据

xarray as xrwrffile = "/home/mw/input/typhoon9537/wrfout_d01_2019-08-08_20_00_00"ds_wrf = xr.open_dataset(wrffile)# 提取降雨量temp = ds_wrf.Ttemp[0,0].plot()

<matplotlib.collections.QuadMesh at 0x7f0c700cfd50>

#meb.tool.plot_tools.contourf_2d_grid(temp)

可试试以上代码，结论是无法直接可视化，需要通过转换投影方式
转换投影方式在workshop有提到，利用xesmf可以将两种不同网格的数据再栅格化到同一网格

2.1.2 网格生成

import xesmf as xeimport xarray as xrimport meteva.base as mebimport numpy as npXLAT = 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





#设计输出网格ds_out = xr.Dataset(                    {                     'lat': (['lat'], np.arange(22, 30+1, 0.1)),                     'lon': (['lon'], np.arange(115, 130+1, 0.1)),                    }                   )

2.1.3 读取era5数据

#meteva读取nc数据并切割lekima = meb.read_griddata_from_nc('/home/mw/input/ERA5_Lekima4742/ERA5_Lekima.nc',value_name='t')lekima.shapelekima =lekima.isel(level=6,time=20)lekima=lekima.sel(lon=slice(116, 131), lat=slice(22, 31))-273.15

2.1.4 重插值

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

Overwrite existing file: patch_437x447_90x160.nc 
 You can set reuse_weights=True to save computing time.

In [29]:

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

<matplotlib.collections.QuadMesh at 0x7f0c70319150>

2.1.5 将插值后的数据放入meteva数组

grid0 = meb.grid([115,130.9,0.1],[22,30.9,0.1],gtime=["2019080820"],dtime_list = [0],level_list = [0],member_list = ["temp"])print(grid0)tnewnew= meb.grid_data(grid0,t_new.data)   #根据网格信息和numpy数组生成网格数print(tnewnew)

members:['temp']
levels:[0]
gtime:['20190808200000', '20190808200000', '1h']
dtimes:[0]
glon:[115, 130.9, 0.1]
glat:[22, 30.9, 0.1]

<xarray.DataArray 'data0' (member: 1, level: 1, time: 1, dtime: 1, lat: 90, lon: 160)>
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 'temp'
  * level    (level) int64 0
  * time     (time) datetime64[ns] 2019-08-08T20:00:00
  * dtime    (dtime) int64 0
  * lat      (lat) float64 22.0 22.1 22.2 22.3 22.4 ... 30.5 30.6 30.7 30.8 30.9
  * lon      (lon) float64 115.0 115.1 115.2 115.3 ... 130.6 130.7 130.8 130.9

2.1.6 完成数据组装后进行绘图

map_extend = [117,129,22,30]axs = meb.creat_axs(2,map_extend,ncol = 2,sup_title = "wrf n era5",                    add_index = ["a","b"],width=8,sup_fontsize = 8)image = meb.add_contourf(axs[0],tnewnew, cmap = "rainbow", clevs = np.arange(0,40,4),add_colorbar=True)image = meb.add_contourf(axs[1],lekima, cmap = "rainbow",clevs = np.arange(0,40,4), add_colorbar=True

)

2.2 nc网格数据插值站点三步走：meteva读取era5 nc格式 地表2m温度数据并插值到气象站点

2.2.1 第一步 读取数据

t2new = meb.read_griddata_from_nc('/home/mw/input/meteva2260/air2m20210905.nc',time='2021-09-05T01:00:00.000000000')

2.2.2 第二步 插值

#格点插到站点 双线性插值

sta = meb.interp_gs_linear(t2new,station)sta

Out[34]:

26544 rows × 7 columns

2.2.3 第三步 画图

meb.tool.plot_tools.scatter_sta(sta[sta['time'] == '2021-09-05 00:00:00'])

结论

1. meteva虽函数简便，但更多图像细化还需matplotlib实现，meteva胜在方便

2. 对于兰伯特投影的wrfout数据无法直接转换，未来可期

3. 自然，一个库无法胜任所有需要。开源库众多正是python的魅力所在。

完整文件与代码后台回复'meteva'，探索不同的库，用不同的、更简洁的方法实现需求，这一过程你会收获良多。