业务刚需！基于meteva的华南降水mode检验

感谢通化的崔忠强大佬来稿，meteva之前一直是用来作图，其实它的本职是搞检验，看看下面看看气象局的业务人员如何使用meteva

MODE方法全称为基于对象诊断评估方法。该方法通过匹配预报与格点实况中的对象，计算出各对象的各类属性，设定不同属性的权重系数，运用模糊逻辑算法计算预报性能的总收益函数从而判断预报的整体表现。可以有效防止传统检验方法在降水落区差异时引起的空报漏报双重惩罚问题。

计算流程如下：

1. 观测场和预报场平滑、雨量阈值提取、选取关注的降雨面积最大最小值，提取目标。

2. 通过目标质心距离、目标面积、最小边界距离等属性，将相似的观测场和预报场目标进行匹配。

3. 观测场和预报场目标存在多对多的情况，将多对多的目标进行合并。

4. 计算目标的多种属性，如凸包点坐标、面积、长短轴长度、方位角等。

5. 通过观测场和预报场的目标属性，计算观测和预报的相似度矩阵。通常认为目标相似度达到0.7表示该模式对雨带的预报较为准确。

本文仅仅是从代码层面简单介绍一下格点实况对格点预报检验的数据处理流程，介绍GFS和CMPAS格点的处理以匹配meteva的检验函数。

具体原理与方法参见meteva官方文档相关页面：https://www.showdoc.com.cn/meteva/6457147666090170

import meteva.base as meb      # 该模块用于IO和基础计算
import meteva.method as mem    # 该模块用于检验的基础算法
import pandas as pd
import numpy as np

#本示例主要对华南地区4月18日的强对流天气进行空间检验。检验数据为GFS模式24小时时效的24小时降水。

#首先初始化检验范围和分辨率。我这边简单采用了metdig给出的华南经纬度范围

grid0 = meb.grid([95,126,0.05],[16,35,0.05])

#首先读取GFS在4月17日20时起报的降水数据
rain24_gfs=meb.read_griddata_from_grib(r".\data\GRAPESGFS_TPE_1_2024041712_NEHE_1_2.grib2",grid=grid0,dtime_dim='step')
rain24_gfs

#这里先看看降水数据是怎么存的。
#给某一点数据做一个随时间演变的切片，输出一下数据
rain24_gfs[0,0,0,:,0,0]

#我们发现这个数据啊，是随预报时效递增的。中间参杂着大量相同的数据。所以我们判断，降水数据是累计的。如果我们要取48小时时效的24小时降水，那么就是48时效的数据减去24时效的数据。

#但是作者偷懒了，想检验的是24小时时效的24小时累计降水，所以直接取dtime在24小时的数据就好了！

rain24_gfs = meb.between_dtime_range(rain24_gfs,start_dtime=24,end_dtime=24)

#为了让数据有辨识度，我们还需要将数据的member属性改为模式名，这样方便程序后续出图的标题处理。

meb.set_griddata_coords(rain24_gfs,member_list = ["CMA-GFS"])

#注意，set_griddata_coords函数是原地变换，也就是直接对rain24_gfs变换，无需接受返回值。

rain24_gfs

#最后把time转换为北京时间，这时候预报场的数据准备完毕啦！

rain24_gfs['time'] = rain24_gfs.time+pd.Timedelta(hours=8)
rain24_gfs

#因为是整点的24小时降水检验，实况场我们采用cmpas的三源融合格点降水实况数据。同样使用meteva去读文件

rain24_cmpas=meb.read_griddata_from_grib(r".\data\Z_SURF_C_BABJ_20240418201539_P_CMPA_FRT_CHN_0P05_DAY-PRE-2024041820.GRB2",grid=grid0)
rain24_cmpas

#这里发现dtime用了时间戳作为数值，我们需要将dtime设置为0.

rain24_cmpas['dtime'] = np.array([0])
#这里同样优化一下member显示
meb.set_griddata_coords(rain24_cmpas,member_list = ["多源融合实况"])
rain24_cmpas

"""

最后我们调用meteva的mode检验主函数。函数包含了目标匹配、相似度计算等mode全流程检验功能。第一参数传入实况场，第二个传入预报场。

smooth对降水（尤其是实况场）做平滑的强度。threshold是识别降水的阈值。这里以大雨以上量级为例，使用25作为识别阈值。

minsize代表了识别出来的目标的最小面积，用来过滤一些叫小尺度的目标。

match_method代表了不同的匹配方法。默认采用centmatch方法。cmap代表了色标。

"""

mem.mode.operate(rain24_cmpas,rain24_gfs,smooth = 5,threshold = 25,minsize = 50,
                    match_method = mem.mode.centmatch, cmap = "rain_24h",
                               save_dir = r"C:\Users\admin\OneDrive\code\tougao",show = True)

"""

{'time': '20240417200000', 'dtime': 24, 'match_count': 2, 'feature_table': {'contingency_table_yesorno': {'Hits': 2, 'Misses': 4, 'False alarms': 1, 'Correct negatives': 93}, 'score': {'ets': 0.2668621700879765, 'pod': 0.3333333333333333, 'pofd': 0.010638297872340425, 'far': 0.3333333333333333, 'hss': 0.42129629629629634}}, 'interester': array([0.89151892, 0.67896435]), 'label_list_matched': [1, 2], 'unmatched': {'ob': [3, 4, 5, 6], 'fo': [3]}, 1: {'feature_axis': {'ob': {'OrientationAngle': {'MajorAxis': 164.44368155700434}, 'window': {'x0': 111.78813952644154, 'y0': 23.001444250045907, 'x1': 103.40841122575851, 'y1': 25.334221273361337}, 'lengths': {'MajorAxis': 7.197021484375, 'MinorAxis': 3.1000823974609375}}, 'fo': {'OrientationAngle': {'MajorAxis': 171.134811755046}, 'window': {'x0': 112.46727240580368, 'y0': 24.358225325285545, 'x1': 104.05592240018089,

...

'feature_props': {'ob': {'centroid': {'x': 105.35, 'y': 26.371785714285714}, 'area': 0.35000000000000003, 'intensity': array([ 9.09829998, 19.89566984, 21.34662037, 27.22632456, 35.09889984, 41.16512394, 47.65501862, 49.94205017, 54.17490005])}}}, 'is_summary': True}

Output is truncated. View as a scrollable element or open in a text editor. Adjust cell output settings...

'''由于我们设置show是True，控制台输出了一坨东西。

其中第一个是匹配目标的空间分布图

第二个图是检验概要表

这时候可以在指定目录下看到如下几个文件夹：

png是图片输出，json是json格式的文本输出，table是匹配结果的表格输出（图片格式）

'''

#完整代码：
import meteva.base as meb      # 该模块用于IO和基础计算
import meteva.method as mem    # 该模块基础了检验的基础算法
import pandas as pd
import numpy as np


grid0 = meb.grid([95,126,0.05],[16,35,0.05])
rain24_gfs=meb.read_griddata_from_grib(r".\data\GRAPESGFS_TPE_1_2024041712_NEHE_1_2.grib2",grid=grid0,dtime_dim='step')
rain24_gfs = meb.between_dtime_range(rain24_gfs,start_dtime=24,end_dtime=24)
meb.set_griddata_coords(rain24_gfs,member_list = ["CMA-GFS"])
rain24_gfs['time'] = rain24_gfs.time+pd.Timedelta(hours=8)
rain24_cmpas=meb.read_griddata_from_grib(r".\data\Z_SURF_C_BABJ_20240418201539_P_CMPA_FRT_CHN_0P05_DAY-PRE-2024041820.GRB2",grid=grid0)
rain24_cmpas['dtime'] = np.array([0])
meb.set_griddata_coords(rain24_cmpas,member_list = ["多源融合实况"])
mem.mode.operate(rain24_cmpas,rain24_gfs,smooth = 5,threshold = 25,minsize = 50,
                    match_method = mem.mode.centmatch, cmap = "rain_24h",
                               save_dir = r"C:\Users\admin\OneDrive\code\tougao",show = True)