利用matplotlib和ERA5数据绘制时间-高度气象综合廓线图
- 1. 效果图
- 2. 总体思路
- 3. 读取数据
- 4. 图形绘制
- 5. 代码完整版
1. 效果图
2. 总体思路
- 气象预报业务中，有种常用的综合廓线图，其本质上是单个站点时间-高度的等高线或者填色图，其中时间是从右到左来看。所以准备好(time, level)的二维数据，然后依次叠加线条和填色就可以，思路很简单，但是绘图中涉及到了很多细节问题，也是琢磨了一阵子，怕以后忘了，记录下学习过程。
# %%
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib.pyplot import MultipleLocator #导入此类，设置坐标轴间隔
import pandas as pd
from scipy.ndimage import gaussian_filter1d
from matplotlib.font_manager import FontProperties
mpl.rcParams["font.size"] = 13
font_cn = FontProperties(fname="C:/Users/xxx/anaconda3/Lib/site-packages/matplotlib/mpl-data/fonts/chinese_ttf/方正楷体简体.ttf", size=14)
- 导入的库主要有三类：
（1）读取数据的xarray和进行简单计算的numpy，pandas.
（2）绘图的matplotlib
（3）平滑曲线的scipy
就是这么简单！
3. 读取数据
- ERA5：这次的例子中使用的是ERA5的再分析数据，首先利用xr.open_dataset读取数据，利用ds.sel()选择需要的时间范围，最后选取需要的变量。图中主要用到的是气温t, 相对湿度r，水平风u,v和垂直速度w。
# read data
data_path = "D:/python/CSDN/data/ERA5-pressure_20230301-0304.nc"
time_range = pd.date_range('2023-02-01', periods=24*3+1, freq='H')
ds = xr.open_dataset(data_path).sel(time=time_range)
rh, temp, u, v, w = ds['r'], ds['t']-273.15, ds['u'], ds['v'], ds['w'] # (time,level, lat, lon)
- 挑选站点和level，原始数据中“time”维度在“level”维度之前，为了绘制综合廓线图，需要用transpose对二者维度的顺序做调换。其次，level是从100-1000hPa, time也是从小到大，为了满足图片中高度从低到高、时间从右到左进行读取的习惯，需要进行reverse的处理。
- 另外一个重要的点是如何从UTC转为北京时。首先用pd.to_datetime(time, utc=True)对time设置为utc, 然后利用.tz_convert(‘Asia/Shanghai’)转换为北京时，最后用.strftime(‘%d%H’).tolist()设置为日+时的显示格式，并设置为列表。琢磨了半天，这种方法最方便！
# %%
# 挑选站点，转换time, level坐标顺序，并调整level（1000hPa->100hPa),time排序（倒序）
lat, lon = 39.8, 116.47
level = [100,150,200,300,400,500,600,700,800,850,900,925,1000]
rh_bj= rh.sel(latitude = lat, longitude = lon, level = level, method = 'nearest').transpose("level","time")[::-1, ::-3].values
temp_bj= temp.sel(latitude = lat, longitude = lon, level = level, method = 'nearest').transpose("level","time")[::-1, ::-3].values
u_bj= u.sel(latitude = lat, longitude = lon, level = level, method = 'nearest').transpose("level","time")[::-1, ::-3].values
v_bj= v.sel(latitude = lat, longitude = lon, level = level, method = 'nearest').transpose("level","time")[::-1, ::-3].values
w_bj= w.sel(latitude = lat, longitude = lon, level = level, method = 'nearest').transpose("level","time")[::-1, ::-3].values
time = rh['time'][::-3]
time_str = pd.to_datetime(time, utc=True).tz_convert('Asia/Shanghai').strftime('%d%H').tolist() #datetime64转为string，和utc+8
- 数据平滑这一步，看自己的需求。如果觉得直接绘制出来的线条还行，这一步可以忽略。如果觉得线条太生硬，可以利用gaussian_filter1d进行平滑，其中sigma越大，表明线条越平滑，相应的也更加失真，自己根据需求可多加尝试。
# 平滑等高线数据
rh_bj = gaussian_filter1d(rh_bj, sigma=0.5) #sigma数值越大 越平滑
w_bj = gaussian_filter1d(w_bj, sigma=1)
temp_bj = gaussian_filter1d(temp_bj, sigma=0.5)
以上是利用ERA5的原始数据绘制综合廓线图的读取数据部分，如果自己有其他模式或者其他格式的数据，只要处理成(level, time)的二维数据，其中level和time坐标都是从大到小排列就好。
4. 图形绘制
- RH填色图：建立画布和坐标后，首先可以利用ax.contourf绘制RH填色图。其中，利用levels指定显示的色标及间隔，extend=’both’添加后colorbar的最大最小值之外还有填色，这种情况下extendrect默认为False，colorbar的两头是三角形，如果不需要三角，可手动设置为True。fig.colorbar中shrink表示colorbar缩短的比例，pad是colorbar离图片的距离，extendfrac=’auto’表示colorbar两头三角形的大小自动和colorbar里每个间隔的长短一致，如果不设置这个，默认的三角形很小，不美观。填色图中唯一遗憾的是cmap没有找到特别合适的，后期还需要自己单独设置colormap。
#建立画布
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(1,1,1)
#添加RH填色图
rh_plot = ax.contourf(rh_bj, cmap = 'Greens',levels=[60,70,80,90,95], extend='both', alpha=0.7)
fig.colorbar(rh_plot, ax=ax, shrink=0.6, pad=0.02,extendfrac='auto') #添加label bar
- temp,w等高线图：等高线的绘制利用ax.contour，比较简单。默认状态下如果数值大于0，则线条为实线，小于0则为虚线，与综合廓线图中的要求一致，所以就不需要单独设置了。ax.clabel用来在等高线中添加label数值，也比较简单。
#添加temp, w等高线
temp_plot = ax.contour(temp_bj, colors='r') # Negative contours default to dashed.
w_plot = ax.contour(w_bj, colors='k',alpha=0.7)
ax.clabel(temp_plot, inline=True, fontsize=10) #添加等高线label
ax.clabel(w_plot, inline=True, fontsize=10)
- 风羽风向标：ax.barbs用来绘制风羽，其中最重要的是需要设置barb_increments。默认的设置里，半个横杆为5，一个横杆为10，一个旗子为50。但是在中国的风羽使用中，半个横杆为2 m/s， 一个为4，旗子为20，所以需要利用barb_increments = dict(half=2, full=4, flag=20)字典的形式进行设置。
- 这里需要注意的是，中国的风羽使用中，flag为空心则为20，实心则为50，但是ax.barbs中旗子要么是空心要么实心，没有办法直接同时设置两种。如果需要设置2种，需要修改源代码。
#画风向标
barb_increments = dict(half=2, full=4, flag=20)
ax.barbs(u_bj, v_bj, color = 'b',length= 7, alpha=0.8, barb_increments=barb_increments)
- 添加横纵坐标对应的时间和高度：y_pos和y_labels分别对应选择显示出来的y轴位置和对应的label，然后利用ax.set_yticks设置就好。x轴的设置同理。此外，ax.xaxis.set_major_locator和ax.xaxis.set_min_locator用来设置主刻度和次刻度的间隔。
#添加横纵坐标对应的时间和高度
#y轴
y_pos = [0,1,3,5,7,10,12]
y_labels = [1000,925,850,700,500,200,100]
ax.set_yticks(y_pos, labels=y_labels) #显示特定xticklabels
ax.set_ylim(-0.5,len(level)+0.1) #上下留出一点空白
#x轴
x_pos = np.arange(len(time_str))
ax.set_xticks(x_pos, labels=time_str)
ax.xaxis.set_major_locator(MultipleLocator(4)) #设置x轴主刻度显示间隔
ax.xaxis.set_minor_locator(MultipleLocator(2)) #设置x轴次刻度显示间隔
- 添加标题、xlabel等图片说明信息：这里的一个小细节是如何使用中文字体。在开头中设置了 font_cn = FontProperties(fname=“C:/Users/xxx/anaconda3/Lib/site-packages/matplotlib/mpl-data/fonts/chinese_ttf/方正楷体简体.ttf”, size=14) ，这里直接利用fontproperties使用就好。想使用什么样的中文字体，下载到指定位置，如上述这样使用即可。
#添加标题、xlabel等图片说明信息
ax.text(23,-2,'Feb',color='r')
ax.text(0,-2,'北京时',fontproperties = font_cn)
ax.text(0,13.3,'Lat:'+str(lat)+' Lon:'+str(lon),size=15,color='r')
ax.text(21,13.3,'000-072',size=15,color='b')
ax.text(10,13.3,'2023020100',size=15,color='b')
ax.text(7,14,'单点3小时间隔综合廓线图', fontproperties = font_cn, size=22,color='b')
ax.set_ylabel("hPa")
大功告成！
5. 代码完整版
# %%
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib.pyplot import MultipleLocator #导入此类，设置坐标轴间隔
import pandas as pd
from scipy.ndimage import gaussian_filter1d
from matplotlib.font_manager import FontProperties
mpl.rcParams["font.size"] = 13
font_cn = FontProperties(fname="C:/Users/xxx/anaconda3/Lib/site-packages/matplotlib/mpl-data/fonts/chinese_ttf/方正楷体简体.ttf", size=14)
# %%
# read data
data_path = "D:/python/CSDN/data/ERA5-pressure_20230301-0304.nc"
time_range = pd.date_range('2023-02-01', periods=24*3+1, freq='H')
ds = xr.open_dataset(data_path).sel(time=time_range)
rh, temp, u, v, w = ds['r'], ds['t']-273.15, ds['u'], ds['v'], ds['w'] # (time,level, lat, lon)
# %%
# 挑选站点，转换time, level坐标顺序，并调整level（1000hPa->100hPa),time排序（倒序）
lat, lon = 39.8, 116.47
level = [100,150,200,300,400,500,600,700,800,850,900,925,1000]
rh_bj= rh.sel(latitude = lat, longitude = lon, level = level, method = 'nearest').transpose("level","time")[::-1, ::-3].values
temp_bj= temp.sel(latitude = lat, longitude = lon, level = level, method = 'nearest').transpose("level","time")[::-1, ::-3].values
u_bj= u.sel(latitude = lat, longitude = lon, level = level, method = 'nearest').transpose("level","time")[::-1, ::-3].values
v_bj= v.sel(latitude = lat, longitude = lon, level = level, method = 'nearest').transpose("level","time")[::-1, ::-3].values
w_bj= w.sel(latitude = lat, longitude = lon, level = level, method = 'nearest').transpose("level","time")[::-1, ::-3].values
time = rh['time'][::-3]
time_str = pd.to_datetime(time, utc=True).tz_convert('Asia/Shanghai').strftime('%d%H').tolist() #datetime64转为string，和utc+8
# 平滑等高线数据
rh_bj = gaussian_filter1d(rh_bj, sigma=0.5) #sigma数值越大 越平滑
w_bj = gaussian_filter1d(w_bj, sigma=1)
temp_bj = gaussian_filter1d(temp_bj, sigma=0.5)
# %%
#建立画布
fig = plt.figure(figsize=(12, 6))
ax = fig.add_subplot(1,1,1)
#添加RH填色图
rh_plot = ax.contourf(rh_bj, cmap = 'Greens',levels=[60,70,80,90,95], extend='both', alpha=0.7)
fig.colorbar(rh_plot, ax=ax, shrink=0.6, pad=0.02,extendfrac='auto') #添加label bar
#添加temp, w等高线
temp_plot = ax.contour(temp_bj, colors='r') # Negative contours default to dashed.
w_plot = ax.contour(w_bj, colors='k',alpha=0.7)
ax.clabel(temp_plot, inline=True, fontsize=10) #添加等高线label
ax.clabel(w_plot, inline=True, fontsize=10)
#画风向标
barb_increments = dict(half=2, full=4, flag=20)
ax.barbs(u_bj, v_bj, color = 'b',length= 7, alpha=0.8, barb_increments=barb_increments)
#添加横纵坐标对应的时间和高度
#y轴
y_pos = [0,1,3,5,7,10,12]
y_labels = [1000,925,850,700,500,200,100]
ax.set_yticks(y_pos, labels=y_labels) #显示特定xticklabels
ax.set_ylim(-0.5,len(level)+0.1) #上下留出一点空白
#x轴
x_pos = np.arange(len(time_str))
ax.set_xticks(x_pos, labels=time_str)
ax.xaxis.set_major_locator(MultipleLocator(4)) #设置x轴主刻度显示间隔
ax.xaxis.set_minor_locator(MultipleLocator(2)) #设置x轴次刻度显示间隔
#添加标题、xlabel等图片说明信息
ax.text(23,-2,'Feb',color='r')
ax.text(0,-2,'北京时',fontproperties = font_cn)
ax.text(0,13.3,'Lat:'+str(lat)+' Lon:'+str(lon),size=15,color='r')
ax.text(21,13.3,'000-072',size=15,color='b')
ax.text(10,13.3,'2023020100',size=15,color='b')
ax.text(7,14,'单点3小时间隔综合廓线图', fontproperties = font_cn, size=22,color='b')
ax.set_ylabel("hPa")
plt.savefig("综合廓线图.jpg")
文章出处登录后可见！
已经登录？立即刷新
