使用python封装绘图函数循环绘制ERA5风场资料的空间分布图
通常，在处理气象海洋资料时，经常会绘制风场的空间分布图进行简单分析，而常常需要连续绘制多天，并将多张子图绘制到同一个图片中，因此这就需要用到循环绘图。
- 同时考虑到下载的ERA5风场资料的经度排列顺序是-180~180°，这里也简单进行了经度转换，将其转换为0 ~ 360的排列顺序。
- 根据每个子图的数据，将选取的时间也在循环中加上
- 考虑到绘制全球的处理时间较久，这里自由选取任意经纬度进行绘制
- 由于每张子图的填色范围是固定的，所以统一绘制colorbar，至于图片最下端
- 箭头为风场风量，填色为风速大小
下面附上代码：
# -*- coding: utf-8 -*- """ Created on %(date)s @author: %(jixianpu)s Email : 211311040008@hhu.edu.cn introduction : keep learning althongh walk slowly """ from matplotlib.ticker import AutoMinorLocator, MultipleLocator, FuncFormatter import numpy as np import pandas as pd import cartopy.crs as ccrs import cartopy.feature as cfeature from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter import matplotlib.pyplot as plt import xarray as xr import glob from scipy import signal path=r'G:/ERA5_200406/' str1='uwnd.200406.nc' str2='vwnd.200406.nc' # b, a = signal.butter(3, [(2/10)/(1/0.25),(2/3)/(1/0.25)], 'bandpass') ######=================== def make_map(ax, title): # set_extent set crs ax.set_extent(box, crs=ccrs.PlateCarree()) land = cfeature.NaturalEarthFeature('physical', 'land', scale, edgecolor='grey', facecolor='grey' ,zorder=2 ) ax.add_feature(land) # set land color ax.coastlines(scale) # set coastline resolution # set coordinate axis ax.set_xticks(np.arange(box[0],box[1]+xstep, xstep),crs=ccrs.PlateCarree()) ax.set_yticks(np.arange(box[2], box[3]+ystep, ystep),crs=ccrs.PlateCarree()) ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label =False))#经度0不加标识 ax.yaxis.set_major_formatter(LatitudeFormatter()) # plt.tick_params(labelsize=25) ax.set_title(title, fontsize=20, loc='left',pad=12) ax.yaxis.set_minor_locator(AutoMinorLocator(5)) ax.xaxis.set_minor_locator(AutoMinorLocator(10)) ax.tick_params(which='minor', direction='out', length=4,width=0.59, right=True, top=True) ax.tick_params(which='major', direction='out', length=8, width=0.99, pad=3, labelsize=10, bottom=True, left=True, right=True, top=True) return ax # # prepare box = [100, 180, -30, 0] scale = '50m' xstep, ystep = 20, 10 # cmap=plt.get_cmap('RdYlBu_r')#'RdYlBu_r' # =======================data============================================ day=(np.arange(19,31)) da = xr.open_dataset(path+str1).sortby("latitude", ascending= True) da2=xr.open_dataset(path+str2).sortby("latitude", ascending= True) #####################reverse################################################ lon_name = 'longitude' # whatever name is in the data da['longitude_adjusted'] = xr.where(da[lon_name] < 0, da[lon_name]%360,\ da[lon_name]) da = ( da .swap_dims({lon_name: 'longitude_adjusted'}) .sel(**{'longitude_adjusted': sorted(da.longitude_adjusted)}) .drop(lon_name)) da = da.rename({'longitude_adjusted': lon_name}) #####################reverse################################################ lon_name = 'longitude' # whatever name is in the data da2['longitude_adjusted'] = xr.where(da2[lon_name] < 0, da2[lon_name]%360,\ da2[lon_name]) da2 = ( da2 .swap_dims({lon_name: 'longitude_adjusted'}) .sel(**{'longitude_adjusted': sorted(da2.longitude_adjusted)}) .drop(lon_name)) da2 = da2.rename({'longitude_adjusted': lon_name}) #####################reverse################################################ fig=plt.figure(figsize=(14,14)) levels=np.arange(0,21,1) # wspace 调整水平的 plt.subplots_adjust(hspace=0.1) plt.tight_layout() count=0 for i in day: count+=1 print(count) tim_s='2004-06-'+str(i)+'-00' # longitude=slice(0,181) u=da.u.sel(level=slice(850,850),latitude=slice(-30,0),longitude=slice(100,180),time=slice(tim_s,tim_s))[0][0] v=da2.v.sel(level=slice(850,850),latitude=slice(-30,0),longitude=slice(100,180),time=slice(tim_s,tim_s))[0][0] w=np.sqrt(u*u+v*v) lon=u.longitude.data lat=u.latitude.data x,y=np.meshgrid(lon,lat) t=u.time step=10 # =======================plot============================================ ax=fig.add_subplot(4,3,count,projection=ccrs.PlateCarree(central_longitude=180)) make_map(ax,pd.to_datetime(t.data).strftime('%Y_%m_%d_%H:00')) ax.quiver(x[::step,::step],y[::step,::step],u.data[::step,::step],v.data[::step,::step],pivot='mid',\ width=0.003,scale=200,headlength=4,headwidth=4, transform=ccrs.PlateCarree(),color='k',angles='xy' ,zorder=2) # ax.set_aspect('auto') plot=ax.contourf(x,y,w,transform=ccrs.PlateCarree(),extend='both',levels=levels ,zorder=1) # cb=fig.colorbar(plot,ax=ax,shrink=0.8,pad=0.05,aspect=15) # cb.ax.tick_params(labelsize=20) # cb.ax.set_title('$°C$',fontsize=20,loc='right') # if count==6: # break ax3=fig.add_axes([0.25,0.1,0.5,0.015]) # 0.25控制距离左边的距离，0.01控制距离下面的距离，最后两位控制color的长度和厚度 cb=fig.colorbar(plot,cax=ax3,shrink=0.9,pad=0.04,ticks=[0,5,10,15,20],aspect=15,orientation='horizontal') cb.ax.tick_params(labelsize=10) cb.ax.set_title('$m/s$',fontsize=15) plt.show() # fig.savefig('wind_850hpa-2004-06-19_30.png',format='png',dpi=100)
绘制结果如下图所示：
