统计数值预报产品生成时间：单个时效

meteorology/devops
目录

之前的文章《统计数值天气预报模式产品生成的典型时间》中介绍如何根据GRIB2文件创建时间统计产品生成的典型时间，并在《Bokeh教程：添加标注》中介绍如何绘制单个时效的统计图形。

本文在第二篇文章的基础上，再一次介绍如何绘制单个时效产品文件的生成情况，并增加额外的统计图形。

声明：本文仅代表作者个人观点，所用数据无法代表真实情况，严禁转载。关于模式系统的相关信息，请以官方发布的信息及经过同行评议的论文为准。

获取数据

载入需要使用到的库

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from tqdm.notebook import trange, tqdm

from nwpc_data.data_finder import find_local_file

生成时间列表，统计截止 2020 年 4 月 29 日之前 60 天的数据

end_date = pd.to_datetime("2020-04-29 00:00:00")
start_date = end_date - pd.Timedelta(days=60)
date_range = pd.date_range(start_date, end_date, freq="D")

获取文件列表

file_list=[
    find_local_file(
        "MODEL_A/grib2/orig",
        start_time=t,
        forecast_time="240h"
    ) 
    for t
    in date_range    
]

获取所有文件的创建时间

from pathlib import Path

pbar = tqdm(total=len(file_list))
ts = []
for f, d in zip(file_list, date_range):
    ts.append(pd.to_datetime(Path(f).stat().st_mtime_ns) - pd.to_datetime(d.date()))
    pbar.update(1)
pbar.close()

s = pd.Series(ts, name="clock", index=date_range)
df = s.to_frame()
df.head()
clock
2020-02-2905:28:28
2020-03-0105:22:41
2020-03-0205:09:02
2020-03-0304:39:53
2020-03-0405:30:53

切尾均值

计算

与上一篇文章一样，计算切尾均值

from scipy import stats
mean_clock = s.mean().ceil("s")
median_clock = s.median()
std = s.std().ceil("s")
std
Timedelta('0 days 00:15:19')

计算切尾均值（0.1）

count = len(df)
trimmed_s = df.sort_values("clock")[int(count*0.1):int(count*0.9)]
trimmed_mean = trimmed_s.mean().loc["clock"].ceil("s")
trimmed_std = trimmed_s.std().loc["clock"].ceil("s")
trimmed_std
Timedelta('0 days 00:07:15')

将 trimmed_mean + std 作为超时阈值，将 trimmed_mean + trimmed_std 作为警告阈值，为数据打标签

trimmed_upper_value = trimmed_mean + trimmed_std
upper_value = trimmed_mean + std

def get_status(clock):
    if clock > upper_value:
        return "late"
    elif clock > trimmed_upper_value:
        return "warn"
    else:
        return "normal"

df["status"] = df["clock"].map(get_status)
df["start_time"] = df.index.map(lambda x: x.strftime("%m/%d"))
df.head()
clockstatusstart_time
2020-02-2905:28:28late02/29
2020-03-0105:22:41late03/01
2020-03-0205:09:02late03/02
2020-03-0304:39:53normal03/03
2020-03-0405:30:53late03/04

绘图

from bokeh.transform import factor_cmap
from bokeh.palettes import Accent3

from bokeh.models.sources import ColumnDataSource
source = ColumnDataSource(df)

late_source = ColumnDataSource(df[df.status=="late"])

from bokeh.models.annotations import Span
from bokeh.models import BoxAnnotation, ColorBar, LabelSet

data = s
p = figure(
    plot_width=1000, 
    plot_height=600,
    y_axis_type="datetime",
    x_axis_type="datetime",
    title="MODEL A GRIB2 ORIG forecast 240h"
)

l = p.circle(
    "index",
    "clock",
    size=5,
    source=source,
    color=factor_cmap(
        'status', 
        palette=Accent3, 
        factors=['normal', 'warn', 'late'],
    ),
)

p.xaxis.axis_label = "date"
p.yaxis.axis_label = "clock"

upper = BoxAnnotation(
    bottom=upper_value,
    fill_alpha=0.1,
    fill_color="yellow",
)
p.add_layout(upper)

trimmed_upper = Span(
    location=trimmed_upper_value, 
    dimension="width", 
    line_color="orange", 
    line_width=1
)
p.add_layout(trimmed_upper)

trimmed_line = Span(
    location=trimmed_mean, 
    dimension="width", 
    line_color="green", 
    line_width=1,
)
p.add_layout(trimmed_line)

labels = LabelSet(
    x="index",
    y="clock",
    text="start_time",
    level="glyph",
    x_offset=5,
    y_offset=5,
    source=late_source,
    render_mode="canvas",
)
p.add_layout(labels)

show(p)

进一步分析

使用 seaborn 对数据进行进一步分析，比较原始数据和截断数据之间的差异。

import seaborn as sns
sns.set(style="whitegrid")

直方图

原始数据直方图

ax = sns.distplot(df["clock"]/pd.Timedelta(minutes=1))
ax.figure.set_size_inches(10, 5)
ax.set_xlabel("clock [minutes from 00:00]")
plt.show()

截断数据直方图

ax = sns.distplot(trimmed_s["clock"]/pd.Timedelta(minutes=1))
ax.figure.set_size_inches(10, 5)
ax.set_xlabel("clock [minutes from 00:00]")
plt.show()

可以看到，无论是原始数据还是截断数据，文件生成时间都属于长尾数据。

箱线图

原始数据

ax = sns.boxplot(
    x=df["clock"]/pd.Timedelta(minutes=1)
)
ax.figure.set_size_inches(10, 5)
ax.set_xlabel("clock [minutes from 00:00]")
plt.show()

截断数据

ax = sns.boxplot(
    x=trimmed_s["clock"]/pd.Timedelta(minutes=1)
)
ax.figure.set_size_inches(10, 5)
ax.set_xlabel("clock [minutes from 00:00]")
plt.show()

可以看到，截断之后，离群点变少了，但中位数还是偏左。

参考

《统计数值天气预报模式产品生成的典型时间》

《Bokeh教程：添加标注》