code详解 | 用python实现气象局降水相态图的绘制

程序由本公众号交流群友友提供，欢迎扫描文末二维码进群！

PS，崔大佬说，如果关注的气象局同仁多的话，他后续会多多分享气象业务相关的代码哟~ 还不快动动小手来点赞转发啦

用途&效果

本文分享两个程序。

• 第一个程序通过meteva和metdig对ec确定性预报进行处理，实现了降水相态空间分布展示的功能；
• 第二个程序通过metdig对ec集合预报数据进行处理，实现了集合预报对单点的降水相态预报分布展示的功能。

效果如下，代码附后。

效果图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()

往期文章推荐

用IDW方法插值站点数据并绘制组图

metpy绘制锋生与冷锋

Python+AI+气象+模式大合集

添加小编微信进入交流群