python--循环绘制ERA5风场的空间分布图

使用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)

绘制结果如下图所示：