code详解 | 用python实现气象局降水相态图的绘制
本帖最后由 自学气象人 于 2024-3-2 15:58 编辑本文作者：通化气象局崔忠强
用途&效果
本文分享两个程序。
第一个程序通过meteva和metdig对ec确定性预报进行处理，实现了降水相态空间分布展示的功能（内嵌开启micaps选项）；
第二个程序通过metdig对ec集合预报数据进行处理，实现了集合预报对单点的降水相态预报分布展示的功能（内嵌开启micaps选项）。
效果如下，代码附后。
效果图1
效果图2
实现代码
code1
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.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='[%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=, 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
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=
#用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=,
var_name='ptype')
init_time = stda_data['time'].values.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': , 'lat': ,'id':},method='nearest')
#dtime的毫秒转小时
data.dtime=data.dtime/3600000000000
#筛选时效
data=data.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 =
counts.append(count)
# 绘制柱状图
fig, ax = plt.subplots(figsize=(10, 6))
width = 0.3
x = np.arange(len(unique_numbers))
y_bottom=
#读取标准的相态色标
cmap= cm_collected.get_cmap('met/precipitation_type_nws')
colors = cmap.colors
#把无降水的颜色改为灰色
colors=
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<i>)
y_bottom+=counts<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()</i></i></i></i></i></i></i>
文章来源于微信公众号：自学气象人
有没有相关的最佳实践或指南可以参考？
页:
[1]