cyeva 是一个由彩云科技天气团队和社区贡献者共同开发的用于对气象要素确定性预报准确率进行快速评测的 Python 开源工具库。
cyeva 将致力于让气象要素确定性预报准确率的自动化评估变得简单直接，将集成常用的确定性预报准确率评估指标，且内部算法广泛使用了 numpy 的向量运算实现，对于大数据量的计算也具有较高的计算效率。
安装
通过pip安装
$ pip install cyeva
注意：由于本项目目前处于beta阶段，并非稳定版本，有可能在后续的发布版中出现不兼容性修改，因此如果是在生产环境安装，安装时建议指定版本号，例如
$ pip install cyeva==0.1.0b0
通过conda安装
cyeva 在 conda-forge 上也做了分发，想用 conda 安装的用户可以执行：
$ conda install -c conda-forge cyeva -y
使用
cyeva 为气温、风和降水编写了专门的对象用于处理对应要素的相关指标。
气温
对于气温这种连续性变量，我们通常会比较关心它的 RMSE 、 MAE 等指标。在 cyeva 中我们可以参照以下的例子来计算此类指标：
import numpy as np
from cyeva import TemperatureComparison
np.random.seed(0) # 指定随机种子以保证得到的结果都是一致的
obs = np.sin(np.arange(100)) * 20 + np.random.random(100) * 5 * np.random.choice([1, -1]) # sin数组叠加随机数组模拟真实气温
fcst = obs + np.random.random(100) * 5 * np.random.choice([1, -1]) # 限制预报在观测的正负5°C以内，这样的样例出来的效果更好一些
temp = TemperatureComparison(obs, fcst, unit='degC')
print('accuracy ration within 1 degC:', temp.calc_diff_accuracy_ratio(limit=1)) # 1度准确率（偏差在1°C以内）
print('accuracy ration within 2 degC:', temp.calc_diff_accuracy_ratio(limit=2)) # 2度准确率（偏差在2°C以内）
print('rss:', temp.calc_rss()) # 剩余平方和
print('rmse:', temp.calc_rmse()) # 均方根误差
print('mae:', temp.calc_mae()) # 平均绝对误差
print('chi square:', temp.calc_chi_square()) # 卡方(χ2)
降水
降水的一个特点是它不具有连续性，因此在对其进行准确性评估的时候，除了通用的指标以外，通常会引入一些其他的指标（例如 TS 评分）。此外，降水有明显的分级特点，还需要按照不同级别的做相应的区别，在 cyeva 中我们可以参照以下的例子来计算用于评估降水的指标：
import numpy as np
from cyeva import PrecipitationComparison
np.random.seed(0) # 指定随机种子以保证得到的结果都是一致的
obs = np.random.random(int(100)) * 50
fcst = np.random.random(int(100)) * 50
precip = PrecipitationComparison(obs, fcst, unit='mm')
print('rss:', precip.calc_rss()) # 剩余平方和
print('rmse:', precip.calc_rmse()) # 均方根误差
print('mae:', precip.calc_mae()) # 平均绝对误差
print('chi square:', precip.calc_chi_square()) # 卡方(χ2)
print('accuracy ratio:', precip.calc_accuracy_ratio()) # 准确率(0级)
print('binary accuracy ratio:', precip.calc_binary_accuracy_ratio()) # 准确率(二分/晴雨)
print('false alarm ratio:', precip.calc_false_alarm_ratio()) # 空报率
print('miss ratio:', precip.calc_miss_ratio()) # 漏报率
print('accuracy ratio:', precip.calc_accuracy_ratio(kind='3h', lev='3')) # 准确率(3小时间隔3级/大雨)
for inv in ['1h', '3h', '12h', '24h']: # 不同间隔下的准确率
for lev in range(7):
lev_str = str(lev)
levp_str = f'+{lev_str}'
print(f'accuracy ratio({inv}|{lev_str}):', precip.calc_accuracy_ratio(kind=inv, lev=lev_str))
if lev > 0:
print(f'accuracy ratio({inv}|{levp_str}):', precip.calc_accuracy_ratio(kind=inv, lev=levp_str))
print('ts:', precip.calc_ts()) # TS评分（默认为1h晴雨TS）
for inv in ['1h', '3h', '12h', '24h']: # 不同间隔下的分级TS评分
for lev in range(7):
lev_str = str(lev)
levp_str = f'+{lev_str}'
print(f'ts({inv}|{lev_str}):', precip.calc_ts(kind=inv, lev=lev_str))
if lev > 0:
print(f'ts({inv}|{levp_str}):', precip.calc_ts(kind=inv, lev=levp_str))
print('ets:', precip.calc_ets()) # ETS评分（默认为1h晴雨ETS）
for inv in ['1h', '3h', '12h', '24h']: # 不同间隔下的分级ETS评分
for lev in range(7):
lev_str = str(lev)
levp_str = f'+{lev_str}'
print(f'ets({inv}|{lev_str}):', precip.calc_ets(kind=inv, lev=lev_str))
if lev > 0:
print(f'ets({inv}|{levp_str}):', precip.calc_ets(kind=inv, lev=levp_str))
print('bias:', precip.calc_bias_score()) # bias评分（默认为1h晴雨bias）
for inv in ['1h', '3h', '12h', '24h']: # 不同间隔下的分级bias评分
for lev in range(7):
lev_str = str(lev)
levp_str = f'+{lev_str}'
print(f'bias({inv}|{lev_str}):', precip.calc_bias_score(kind=inv, lev=lev_str))
if lev > 0:
print(f'bias({inv}|{levp_str}):', precip.calc_bias_score(kind=inv, lev=levp_str))
风
对于风这种矢量要素，我们需要同时提供速度和方向信息，因此在实例化对象的时候传入的数据数组会和气温、降水不一样，同时也有一些专门针对于风评估的指标，例如风级偏强率偏弱率等，在 cyeva 中我们可以参照以下的例子来计算用于评估风的指标：
import numpy as np
from cyeva import WindComparison
np.random.seed(0)
obs_spd = np.random.random(100) * 10
obs_dir = np.random.random(100) * 360
fct_spd = np.random.random(100) * 10
fct_dir = np.random.random(100) * 360
wind = WindComparison(obs_spd, fct_spd, obs_dir, fct_dir)
print('difference accuracy ratio within 1 m/s:', wind.calc_diff_accuracy_ratio(limit=1)) # 1m/s准确率（风速偏差在1m/s以内）
print('difference accuracy ratio within 2 m/s:', wind.calc_diff_accuracy_ratio(limit=2)) # 2m/s准确率（风速偏差在2m/s以内）
print('wind speed rss:', wind.calc_rss()) # 剩余平方和（默认风速）
print('wind direction rss:', wind.calc_rss(kind='direction')) # 剩余平方和（指定风向）
print('wind speed rmse:', wind.calc_rmse()) # 均方根误差（默认风速）
print('wind direction rmse:', wind.calc_rmse(kind='direction')) # 均方根误差（指定风向）
print('wind speed mae:', wind.calc_mae()) # 平均绝对误差（默认风速）
print('wind direction mae:', wind.calc_mae(kind='direction')) # 平均绝对误差（指定风向）
print('wind speed chi square:', wind.calc_chi_square()) # 卡方(χ2)
print('wind direction chi square:', wind.calc_chi_square(kind='direction')) # 卡方(χ2)（指定风向）
print('wind direction score:', wind.calc_dir_score()) # 风向评分
print('wind speed score:', wind.calc_speed_score()) # 风速评分
print('wind scale accuracy ratio:', wind.calc_wind_scale_accuracy_ratio()) # 风级准确率
print('wind speed accuracy ratio within 2m/s:', wind.calc_speed_accuracy_ratio()) # 风速准确率(默认2m/s偏差以内)
print('wind speed accuracy radio within 3m/s:', wind.calc_speed_accuracy_ratio(limit=3)) # 风速准确率(指定3m/s偏差以内)
print('wind scale stronger ratio:', wind.calc_wind_scale_stronger_ratio()) # 风级偏强率
print('wind scale weaker ratio:', wind.calc_wind_scale_weaker_ratio()) # 风级偏弱率
其他连续性要素
对于像气压、相对湿度等连续性要素来说，其测评指标通常与气温相同。但是为了代码的优雅，我们在代码里会预留一个通用的 Comparison 对象用于对通用要素做测评，以相对湿度为例，代码示例：
import numpy as np
from cyeva import Comparison
np.random.seed(0) # 指定随机种子以保证得到的结果都是一致的
obs = np.random.random(100) * 100
fcst = obs + np.random.random(100) * 30 * np.random.choice([1, -1])
fcst[fcst>100] = 100
fcst[fcst<0] = 0
temp = Comparison(obs, fcst)
print('accuracy ration within 10 %:', temp.calc_diff_accuracy_ratio(limit=10)) # 10百分比准确率（偏差在10%以内）
print('accuracy ration within 20 %:', temp.calc_diff_accuracy_ratio(limit=20)) # 20百分比准确率（偏差在20%以内）
print('rss:', temp.calc_rss()) # 剩余平方和
print('rmse:', temp.calc_rmse()) # 均方根误差
print('mae:', temp.calc_mae()) # 平均绝对误差
print('chi square:', temp.calc_chi_square()) # 卡方(χ2)
一些实用函数
由于在计算测评指标时，不可避免地需要对数据进行一些转换处理，因此 cyeva 中包含了很多实用的要素处理及转换的函数。我们对这些函数都做了性能优化，既能保证其功能实现，又能保证其运行效率。
-
风速定级：根据风速值确定风级，可以省去查风速表的麻烦。
from cyeva import identify_wind_scale
print(identify_wind_scale(5.5)) # 结果为4
-
风向方位转换：根据角度值计算风的方位，支持8方位和16方位。
from cyeva import identify_direction
# 默认8方位
print(identify_direction(0)) # 0°的风向方位，结果为0（正北）
print(identify_direction(45)) # 45°的风向方位，结果为1（东北）
print(identify_direction(90)) # 90°的风向方位，结果为2（正东）
# 以此类推
# 16方位
print(identify_direction(0, dnum=16)) # 0°的风向方位，结果为0（正北）
print(identify_direction(45, dnum=16)) # 45°的风向方位，结果为2（东北）
print(identify_direction(90, dnum=16)) # 90°的风向方位，结果为4（正东）
返回编号与方位的关系如下表：
-
计算角度偏移量：根据输入的两个角度值计算角度最小偏移量
from cyeva import get_least_angle_deflection
print(get_least_angle_deflection(20, 85)) # 20°与85°的偏移量，结果为65
print(get_least_angle_deflection(0, 355)) # 0°与355°的偏移量，结果为5
print(get_least_angle_deflection(0, 355, absolute=False)) # 0°与355°的偏移量（不取绝对值），结果为-5
性能测试
我们基于 Continuous Benchmark 框架对 cyeva 的部分函数做了持续的性能测试。
benchmark网址：https://caiyunapp.github.io/cyeva/performance/
以 TS 评分的计算为例，它的性能表现如下：
| 数据量级 | 每秒执行次数 | 单次执行耗时(s) |
|---|---|---|
| 1E3 | ~3400 | ~0.0003 |
| 1E6 | ~5 | ~0.2 |
| 1E7 | ~0.5 | ~2 |
欢迎社区的代码贡献
本项目已基于 Apache 协议开源，欢迎广大社区开发者为项目提交 PR 或者 Issue。项目地址：https://github.com/caiyunapp/cyeva。
另外，本项目已经编写了一份说明文档（地址：https://cyeva.readthedocs.io/zh-cn/latest/index.html），文档内包含各种指标的算法说明，欢迎大家的查阅，也欢迎对文档进行补充。文档基于 sphinx 架构编写，源代码集成在仓库的 docs 目录内，且与 Readthedocs 做了集成，欢迎提交 PR。
原文始发于微信公众号（Clarmy吱声）：cyeva: 一个通用的气象预报准确率测评工具包
免责声明:文章中涉及的程序(方法)可能带有攻击性，仅供安全研究与教学之用，读者将其信息做其他用途，由读者承担全部法律及连带责任，本站不承担任何法律及连带责任；如有问题可邮件联系(建议使用企业邮箱或有效邮箱,避免邮件被拦截，联系方式见首页)，望知悉。
- 左青龙
- 微信扫一扫
-
- 右白虎
- 微信扫一扫
-
评论