python气象数据可视化学习记录1——基于ERA5数据画风场和海平面气压填色叠加图
1. 写在前面
大学时为了下载某些资源申请CSDN账号,11年过去了,一篇博文也没有,惭愧。受小王同学的启发,应该把平时自己的一些东西和学习经验记录下来,也是一种很好的经验。所以有了这第一篇文章。
学习python想了好几年了,但因为自己受NCL影响太多,一直走不了这个坑。近些年来,python在气象和大气环境领域有了很好的发展,自己也有了危机感,也决定有一些改变。记录博文也是为了督促自己,当然如果有人能从中得到收获,那也是更开心的。
因为自己还是新手,很多东西是看书+自己琢磨+官网+各类公众号和气象家园的综合,如果有任何错误和版权问题,欢迎大家指出!
2. 图片效果
最近自己科研工作上有一个需求,要利用ERA5 0.25°的小时资料画风场和海平面气压叠加的图,包含4个子图,分别代表了不同的时间,效果如下:
- 填色代表海平面气压hPa,箭头为风场
3. 逐步代码解析
3.1导入库
使用到的库主要有
- matplotlib用来做图,包含了画图的pyplot,ticker等函数
- numpy 简单计算
- cartopy 绘地图
- netCDF4 读取NC格式数据
import matplotlib.pyplot as plt###引入库包
import numpy as np
import matplotlib as mpl
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker
import netCDF4 as nc
import matplotlib.pyplot as plt
此外,为了控制图内所有字体大小、轴线粗细、字体等信息一致,使图片更加美化,统一进行如下设置:
mpl.rcParams["font.family"] = 'Arial' #默认字体类型
mpl.rcParams["mathtext.fontset"] = 'cm' #数学文字字体
mpl.rcParams["font.size"] = 12 #字体大小
mpl.rcParams["axes.linewidth"] = 1 #轴线边框粗细(默认的太粗了)
最后,一个python大神告诉我的方法,加入一个魔法函数,可以不使用plt.show()就可以显示图片的方法:
%matplotlib inline
到此,引入库包的步骤就完成啦!
3.2 读取NC格式数据
读取文件,这里用了ERA5的pressure文件来选取850hpa上的uv风场,single文件里则代表了地表的参量。
f1=nc.Dataset('E:/winter-photo-data/sum/data/ERA-data/202010-ERA-single.nc')
f2=nc.Dataset('E:/winter-photo-data/sum/data/ERA-data/202010-ERA-pressure.nc')
读取变量,需要用到msl(海平面气压),lat,lon,time,u,v等信息。
注意:
- 如果不知道变量名,可以print(f1)来看下具体有哪些变量。
- ERA5中lat是从50->30这样从北往南排列的。
- lat,lon 都是一维函数,msl是(时间,lat,lon)三维函数,u850和v850都是(时间,层数,lat,lon)的四位函数,绘图时只需要二维结果,所以后续要对数据进行挑选和加工。
msl=f1.variables["msl"][:]/100 #pa->hpa
lat=f1.variables['latitude'][:]
lon=f1.variables['longitude'][:]
time=f1.variables['time'][:]
u850=f2.variables["u"][:]
v850=f2.variables["v"][:]
3.3 对数据进行加工1
首先,要建立二维的lon2d,和lat2d,后面画填色图的时候经纬度信息要是二维的。
lon2d, lat2d = np.meshgrid(lon, lat)
其次,用布尔索引的方法挑选需要的时间段。
time_all=[(time>=1058760) & (time<=1059096)] #13-27 Oct
time_20=[(time>=1058928) & (time<=1058951)] # 20 Oct
time_25=[(time>=1059048) & (time<=1059071)] # 25 Oct
time_clean=[(time>=1058760) & (time<1058928)]
这样,数据就挑选完成了,接下来将用lon2d,lat2d,msl,u850,v850,还有time等逻辑值来画图啦。
3.4 建立画布和子图,选择投影方式
使用的是最简单的lat,lon的投影,图片大小是10*8,有两行两列4个子图,这样就建立了画布,得到了各个子图的轴ax。
proj = ccrs.PlateCarree(central_longitude=130) #中国为左
fig = plt.figure(figsize=(10,8),dpi=550) # 创建画布
ax = fig.subplots(2, 2, subplot_kw={'projection': proj}) # 创建子图
3.5 (最重要的部分)自定义画图函数
3.5.1 为什么要定义画图函数?
要画多张图,可统一定义函数,后面直接传递参数,调用画图即可。不然要重复多次步骤。
3.5.2 生成二维的风场和海平面气压数据
定义的函数名为plot4,传递的参量只有time_range(对应time_all, time_20等布尔逻辑值,用来挑选所需要的时次),ax轴,labels为图左上角的标签。
msl_all等为时间平均后的二维数组,u_all等为时间平均后包含高度层信息的三维数组,后面再直接选择特定的层就行,这里没有选择。
def plot4(time_range,ax,labels):
msl_all=np.mean(msl[time_range],axis=0)
u_all=np.mean(u850[time_range],axis=0)
v_all=np.mean(v850[time_range],axis=0)
3.5.3 利用cartopy绘制地图底图
这里参考了小白学习cartopy画地图的第一天(中国行政区域图,含南海)的博文,添加了经纬度信息,并对中国地区的边界进行了重新的设定。具体不再说明,上面的文章里写的很详细。如果用cartopy画地图的话可以学习下。需要注意的是选择extent=[100,125,30,50]经纬度信息时,一定要跟你画图时二维数组的经纬度信息一致,不然出的图是错的。
#-----------绘制地图-------------------------------------------
ax.add_feature(cfeature.LAND.with_scale('50m'))####添加陆地######
ax.add_feature(cfeature.COASTLINE.with_scale('50m'),lw=0.3)#####添加海岸线#########
ax.add_feature(cfeature.OCEAN.with_scale('50m'))######添加海洋########
#-----------添加经纬度---------------------------------------
extent=[100,125,30,50]##经纬度范围
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0., color='k', alpha=0.5, linestyle='--')
gl.top_labels = False ##关闭上侧坐标显示
gl.right_labels = False ##关闭右侧坐标显示
gl.xformatter = LONGITUDE_FORMATTER ##坐标刻度转换为经纬度样式
gl.yformatter = LATITUDE_FORMATTER
gl.xlocator = mticker.FixedLocator(np.arange(extent[0], extent[1]+1, 5))
gl.ylocator = mticker.FixedLocator(np.arange(extent[2], extent[3]+1, 5))
gl.xlabel_style={'size':10}
gl.ylabel_style={'size':10}
ax.set_extent(extent) #显示所选择的区域
#----------修改国界,并添加省界-----------------------------
with open('C:/Users/hj/.local/share/cartopy/shapefiles/natural_earth/physical/CN-border-La.dat') as src:
context = src.read()
blocks = [cnt for cnt in context.split('>') if len(cnt) > 0]
borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks]
for line in borders:
ax.plot(line[0::2], line[1::2], '-', color='k',lw=0.3, transform=ccrs.Geodetic())
3.5.4 画填色图和风场图
- 首先levels定义了填色图色标的范围,1004-1032之间,间隔1个给一个颜色。
- ax.contourf绘制了海平面气压的填色图,参数lon2d, lat2d ,msl_all都是二维的数组,colormap选择了’Spectral_r’,这个可以看个人的需求,网上有很多资源。
- ax.quiver是绘制风场的,也需要输入lon2d, lat2d,u,v的信息(此时uv为三维数组,需要挑选自己所需要的层数,变成二维的),scale用来调整画出来的风场的大小,值越大,风的线条越小。width是风线条的粗细,headwidth和headlength是箭头粗细和长度。uv和lat,lon都是每隔5个显示,是为了使风场图不那么密集,这个可以根据自己的需求进行调节。
最重要的是,一定要加上transform=ccrs.PlateCarree(),图片中地图需要转换下,使地图和数据直接做出来的二维图对应上。 - 绘制每个子图的标题,pad代表标题离轴的远近。
返回了填色图和风场图,这样函数就完成了,后面只需要传递参数调用就行。
#-------------------plot---------------------------
levels = np.arange(1004, 1032 + 1, 1)
cb = ax.contourf(lon2d,lat2d,msl_all, levels=levels, cmap='Spectral_r',transform=ccrs.PlateCarree())
cq = ax.quiver(lon2d[::5,::5],lat2d[::5,::5],u_all[16,::5,::5],v_all[16,::5,::5],color='k',
scale=75,zorder=10,width=0.005,headwidth=3,headlength=4.5,transform=ccrs.PlateCarree()) #925hpa->16
ax.set_title(labels,loc="left",fontsize=12,pad=1)
return cb, cq
3.6 调用函数,画图
subplots_adjust调整了图片的大小和子图之间的距离等,使图片更加美观。
后面直接进行4次函数调用,就可以啦。
#------------------------plot-----------------------
plt.subplots_adjust(left=0.15,right=0.85,top=0.8,bottom=0.2,wspace=0.15,hspace=0.2)
cb1,cq1=plot4(time_all,ax[0][0],'(a) Mean')
cb2,cq2=plot4(time_clean,ax[0][1],'(b) Clean days')
cb3,cq3=plot4(time_20,ax[1][0],'(c) 10/20')
cb4,cq4=plot4(time_25,ax[1][1],'(d) 10/25')
3.7 添加colorbar
因为是填色图,所以一定要条件colorbar. 因为是4幅图一个色标,所以ax=(ax[0][0],ax[0][1],ax[1][0],ax[1][1]),方向为垂直放置在cb3旁边,tick显示的范围为1004-1032,每隔4个显示。aspect是colorbar的长宽高比例,shrink是代表要显示多长,是百分比。pad是离图片的距离。
cbar=fig.colorbar(cb3,ax=(ax[0][0],ax[0][1],ax[1][0],ax[1][1]),orientation='vertical',ticks=np.arange(1004, 1032 + 4,4),
aspect=30,shrink=0.9,pad=0.02)
cbar.ax.tick_params(pad=0.01,length=2,width=0.15)
ax.tick_params里的pad的意思是色标上的数字离色标的距离,length和width是色标tick的长断粗细。
3.8 添加图中的text
为了在图中标出高压和低压中心,需要用ax.text的进行添加。ax为需要添加到哪个图中,0.85和0.87则用百分比的方式决定了点的位置,其他不再赘述,比较简单。一定要添加transform=ax[1][0].transAxes,进行坐标转换。
#--------------------add text----------------------------
ax[1][0].text(0.85,0.87,'L',color='red',fontsize=15,weight='bold',va='bottom',ha='center',transform=ax[1][0].transAxes)
ax[1][1].text(0.92,0.12,'H',color='red',fontsize=15,weight='bold',va='bottom',ha='center',transform=ax[1][1].transAxes)
3.9 图片输出
可以输出eps,pdf等矢量图,也可以有Png这样的栅格图。OK,大功告成!
plt.savefig("Figure2.png")
4. 代码完整版
import matplotlib.pyplot as plt###引入库包
import numpy as np
import matplotlib as mpl
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker
import netCDF4 as nc
import matplotlib.pyplot as plt
mpl.rcParams["font.family"] = 'Arial' #默认字体类型
mpl.rcParams["mathtext.fontset"] = 'cm' #数学文字字体
mpl.rcParams["font.size"] = 12
mpl.rcParams["axes.linewidth"] = 1
%matplotlib inline
f1=nc.Dataset('E:/winter-photo-data/sum/data/ERA-data/202010-ERA-single.nc')
f2=nc.Dataset('E:/winter-photo-data/sum/data/ERA-data/202010-ERA-pressure.nc')
msl=f1.variables["msl"][:]/100 #pa->hpa
lat=f1.variables['latitude'][:]
lon=f1.variables['longitude'][:]
time=f1.variables['time'][:]
u850=f2.variables["u"][:]
v850=f2.variables["v"][:]
lon2d, lat2d = np.meshgrid(lon, lat)
#-----------------mean-----------------------------
time_all=[(time>=1058760) & (time<=1059096)] #13-27 Oct
time_20=[(time>=1058928) & (time<=1058951)] # 20 Oct
time_25=[(time>=1059048) & (time<=1059071)] # 25 Oct
time_clean=[(time>=1058760) & (time<1058928)]
proj = ccrs.PlateCarree(central_longitude=130) #中国为左
fig = plt.figure(figsize=(10,8),dpi=550) # 创建画布
ax = fig.subplots(2, 2, subplot_kw={'projection': proj}) # 创建子图
def plot4(time_range,ax,labels):
msl_all=np.mean(msl[time_range],axis=0)
u_all=np.mean(u850[time_range],axis=0)
v_all=np.mean(v850[time_range],axis=0)
#-----------绘制地图-------------------------------------------
ax.add_feature(cfeature.LAND.with_scale('50m'))####添加陆地######
ax.add_feature(cfeature.COASTLINE.with_scale('50m'),lw=0.3)#####添加海岸线#########
ax.add_feature(cfeature.OCEAN.with_scale('50m'))######添加海洋########
#-----------添加经纬度---------------------------------------
extent=[100,125,30,50]##经纬度范围
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0., color='k', alpha=0.5, linestyle='--')
gl.top_labels = False ##关闭上侧坐标显示
gl.right_labels = False ##关闭右侧坐标显示
gl.xformatter = LONGITUDE_FORMATTER ##坐标刻度转换为经纬度样式
gl.yformatter = LATITUDE_FORMATTER
gl.xlocator = mticker.FixedLocator(np.arange(extent[0], extent[1]+1, 5))
gl.ylocator = mticker.FixedLocator(np.arange(extent[2], extent[3]+1, 5))
gl.xlabel_style={'size':10}
gl.ylabel_style={'size':10}
ax.set_extent(extent) #显示所选择的区域
#----------修改国界,并添加省界-----------------------------
with open('C:/Users/hj/.local/share/cartopy/shapefiles/natural_earth/physical/CN-border-La.dat') as src:
context = src.read()
blocks = [cnt for cnt in context.split('>') if len(cnt) > 0]
borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks]
for line in borders:
ax.plot(line[0::2], line[1::2], '-', color='k',lw=0.3, transform=ccrs.Geodetic())
#-------------------plot---------------------------
levels = np.arange(1004, 1032 + 1, 1)
cb = ax.contourf(lon2d,lat2d,msl_all, levels=levels, cmap='Spectral_r',transform=ccrs.PlateCarree())
cq = ax.quiver(lon2d[::5,::5],lat2d[::5,::5],u_all[16,::5,::5],v_all[16,::5,::5],color='k',
scale=75,zorder=10,width=0.005,headwidth=3,headlength=4.5,transform=ccrs.PlateCarree()) #925hpa->16
ax.set_title(labels,loc="left",fontsize=12,pad=1)
return cb, cq
#------------------------plot-----------------------
plt.subplots_adjust(left=0.15,right=0.85,top=0.8,bottom=0.2,wspace=0.15,hspace=0.2)
cb1,cq1=plot4(time_all,ax[0][0],'(a) Mean')
cb2,cq2=plot4(time_clean,ax[0][1],'(b) Clean days')
cb3,cq3=plot4(time_20,ax[1][0],'(c) 10/20')
cb4,cq4=plot4(time_25,ax[1][1],'(d) 10/25')
cbar=fig.colorbar(cb3,ax=(ax[0][0],ax[0][1],ax[1][0],ax[1][1]),orientation='vertical',ticks=np.arange(1004, 1032 + 4,4),
aspect=30,shrink=0.9,pad=0.02)
cbar.ax.tick_params(pad=0.01,length=2,width=0.15)
#--------------------add text----------------------------
ax[1][0].text(0.85,0.87,'L',color='red',fontsize=15,weight='bold',va='bottom',ha='center',transform=ax[1][0].transAxes)
ax[1][1].text(0.92,0.12,'H',color='red',fontsize=15,weight='bold',va='bottom',ha='center',transform=ax[1][1].transAxes)
plt.savefig("Figure2.png")