程序由本公众号交流群友友提供，欢迎扫描文末二维码进群！
PS，崔大佬说，如果关注的气象局同仁多的话，他后续会多多分享气象业务相关的代码哟~ 还不快动动小手来点赞转发啦
用途&效果
本文分享两个程序。
第一个程序通过meteva和metdig对ec确定性预报进行处理，实现了降水相态空间分布展示的功能（内嵌开启micaps选项）； 第二个程序通过metdig对ec集合预报数据进行处理，实现了集合预报对单点的降水相态预报分布展示的功能（内嵌开启micaps选项）。
效果如下，代码附后。
效果图1
效果图2
实现代码
code1
"""
@author: 通化市气象局崔忠强
"""
from datetime import datetime
import metdig
metdig.set_loglevel('debug')
import metdig.graphics.cmap.cm as cm_collected
from metdig.graphics.boxplot_method import *
from metdig.graphics.draw_compose import *
import metdig.graphics.cmap.cm as cm_collected
import matplotlib.pyplot as plt
import meteva.base as meb
from metdig.graphics.pcolormesh_method import *
#利用meteva读取ec的降水相态确定性预报grib2文件
data=meb.read_griddata_from_grib(r"C:\Users\admin\OneDrive\code\meteva_jianyan\era\ECMFC1D_PRTY_1_2024011600_GLB_1_2.grib2",grid=None,dtime_dim='step')
#手动输入起报时间
init_time = data['time'].values[0].tolist()
init_time = datetime.fromtimestamp(init_time/1000000000)
print(init_time)
#输入绘图范围
map_extent = (100,130,30,55)
#输入预报时效
fhour = 48
init_time_str=init_time.strftime('%Y%m%d%H%M')
data = data.sel(dtime=fhour)
#如果有接入micaps服务器的条件，可以调用micaps服务器的数据，这样可以省去上面预处理本地grib的步骤：
#from metdig.io.cassandra import get_model_grids,get_model_grid
#init_time=datetime(2024,1,16,8)
#data = get_model_grid(data_name='ecmwf', var_name='ptype', init_time=init_time, fhour=fhour)
#利用metdig初始化画布
obj = horizontal_compose(title='[ECWMF][%s] 起报 [%d] 时效降水相态'%(init_time_str,fhour), description='[%s] 起报 [%d] 时效降水相态'%(init_time_str,fhour), png_name='降水相态', map_extent=map_extent, is_return_figax=True,figsize=(16, 12))
#使用metdig的pcolormesh_2d绘制
img=pcolormesh_2d(obj.ax, data, xdim='lon', ydim='lat',
add_colorbar=False, cb_pos='bottom', cb_ticks=None, cb_label=None,
levels=[0,1,3,5,6,7,8], cmap= cm_collected.get_cmap('met/precipitation_type_nws'), extend='both',
transform=ccrs.PlateCarree(), alpha=1)
#添加色标
cb = plt.colorbar(img, ax=obj.ax, shrink=0.8, pad=0.05,extend='max')
# 设置色标文字标签
labels = ['无', '雨', '冻雨', '雪', '雨夹雪', '湿雪', '冰粒']
cb.ax.set_yticklabels(labels)
plt.show()
code2
"""
@author: 通化市气象局崔忠强
"""
from datetime import datetime
import metdig
metdig.set_loglevel('debug')
import metdig.graphics.cmap.cm as cm_collected
from metdig.graphics.boxplot_method import *
from metdig.graphics.draw_compose import *
import metdig.graphics.cmap.cm as cm_collected
import matplotlib.pyplot as plt
from metdig.graphics.pcolormesh_method import *
import metdig.utl as mdgstda
import copy
import cfgrib
#初始化要素
lon=113
lat=34
id='郑州'
fhours=[24,30,36,42,48,54,60,66,72,78,84,96]
#用cfgrib读取集合预报
data = cfgrib.open_dataset(r"C:\Users\admin\OneDrive\code\meteva_jianyan\era\ECMFENS_PRTY_1_2024011600_ASI_4_4.grib2")
#转换成metdig的格式，方便调用metdig进行插值
member = data.coords['number'].values
member = list(map(lambda x: 'ecens' + '-' + str(x), member))
stda_data = None
stda_data = mdgstda.xrda_to_gridstda(data['ptype'],
member_dim='number', level_dim='surface', time_dim='time', lat_dim='latitude', lon_dim='longitude',dtime_dim='step',
member=member, level=[0],
var_name='ptype')
init_time = stda_data['time'].values[0].tolist()
init_time = datetime.fromtimestamp(init_time/1000000000)
print(init_time)
init_time_str=init_time.strftime('%Y%m%d%H%M')
#如果有接入micaps服务器的条件，可以调用micaps服务器的数据，这样可以省去上面预处理本地grib的步骤：
#from metdig.io.cassandra import get_model_grids,get_model_grid
#init_time = datetime(2024, 1, 16, 0)
#data = get_model_grids(data_name='ecmwf_ens', var_name='ptype', init_time=init_time, fhours=fhours)
#插值到站点（临近点插值）
data=mdgstda.gridstda_to_stastda(stda_data, points={'lon': [lon], 'lat': [lat],'id':[id]},method='nearest')
#dtime的毫秒转小时
data.dtime=data.dtime/3600000000000
#筛选时效
data=data[data['dtime'].isin(fhours)]
# 创建示例数据
timedelta = pd.to_timedelta(data['dtime'], unit='h')
data['new_time'] = pd.to_datetime(data['time']) + timedelta + pd.to_timedelta(8, unit='h')#北京时间
times = data['new_time'].tolist()
#times的每个元素去掉前面的年月
for i in range(len(times)):
times[i] = times[i].strftime('%d %H:%M')
#删除data的前6列基础数据
data = data.drop(data.columns[:6], axis=1)
#删除data的最后一列time new
data = data.drop(data.columns[-1], axis=1)
#把data转换为整型方便做映射
data = data.astype(int)
#把data中的5改为文字 雪
data = data.replace(5, '雪')
data = data.replace(0, '无')
data = data.replace(1, '雨')
data = data.replace(3, '冻雨')
data = data.replace(6, '湿雪')
data = data.replace(7, '雨夹雪')
data = data.replace(8, '冰粒')
data = data.replace(12, '雨夹雪')
data = data.values.tolist()
#计算每个相态出现次数
unique_numbers = ['无', '雨', '冻雨', '雪', '雨夹雪', '湿雪','冰粒']
counts = []
for i in range(len(times)):
count = [data[i].count(number) for number in unique_numbers]
counts.append(count)
# 绘制柱状图
fig, ax = plt.subplots(figsize=(10, 6))
width = 0.3
x = np.arange(len(unique_numbers))
y_bottom=[0 for i in range(len(times))]
#读取标准的相态色标
cmap= cm_collected.get_cmap('met/precipitation_type_nws')
colors = cmap.colors
#把无降水的颜色改为灰色
colors[0]=[0.8,0.8,0.8]
for i, num in enumerate(unique_numbers):
y_bottom_old=copy.deepcopy(y_bottom)
y_num=[]
for j in range(len(counts)):
y_num.append(counts[j][i])
y_bottom[j]+=counts[j][i]
ax.bar(times, y_num, width=width, align='center', label=unique_numbers[i], bottom=y_bottom_old,color=colors[i])
# 设置图表标题和轴标签
ax.set_title('EC集合预报[%s]起报[%s]降水相态分布'%(init_time_str,id))
ax.set_xlabel('time')
ax.set_ylabel('count')
# 设置x轴刻度标签
ax.set_xticks(times)
ax.set_xticklabels(times)
# 添加图例
ax.legend()
# 显示图表
plt.show()
往期文章推荐
添加小编微信进入交流群
