python绘图--由逐日风场数据计算月平均风场数据并绘制二维填色图

前言

• 如下图所示，现有1月逐日的一天四次的风场资料，现需要计算1月月平均数据并绘制二维空间填色图
• 处理思路如下：
  • 1、通过for循环读取数据
  • 2、对于读取的日数据按照时间维求平均，将一天四次的资料处理为日平均资料
  • 3、读取风速分量，将一个月内所以日平均资料存到一个空矩阵中
  • 4、计算存储的风速矩阵，按照时间维求平均得到月平均数据
  • 5、将月平均数据进行绘图
• 注意点：
• 对于数据名称比较统一，使用字符串进行循环读取
• 考虑到数据大小，可以在读取数据时截取部分区域，节省时间
• 先读取一个数据查看要设置的矩阵的size
  

code:

# -*- coding: utf-8 -*-
"""
Created on %(date)s

@author: %(jixianpu)s

Email : 211311040008@hhu.edu.cn

introduction : keep learning althongh walk slowly
"""
##########导入库###############################################################
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter,LatitudeFormatter
# import glob

##########导入文件、创建空矩阵用来存储日数据################################


str0 = 'D:\mql\wind\\'
str1 = 'CCMP_Wind_Analysis_200001'
str2 = '_V02.0_L3.0_RSS.nc'



u_mean=np.zeros((160,200,31))

u_mean[:] = np.nan

v_mean=np.zeros((160,200,31))

v_mean[:] = np.nan



##########循环读入日数据、将存储的日数据求平均，得到月平均########################
i = 0

for i in range(1,31):
    file=str0+str1+'%02.0f'%(i)+str2
    
    
    da=xr.open_dataset(file).sel(longitude=slice(130,180),latitude=slice(-20,20)) #按照经纬度区域读取数据

    da_m = da.mean('time')  #  对数据的时间那一维求平均，使得一天四次的数据为一天一次
    lon  = da.longitude.data
    lat  = da.latitude.data
    u_mean[:,:,i] =da_m.uwnd  #  读取u 分量
    v_mean[:,:,i] =da_m.uwnd  #  读取v 分量
    
    x,y=np.meshgrid(lon.data,lat.data)
    i=i+1
#################计算月平均数据，按照时间纬度计算################################
u_monthlymean =np.nanmean(u_mean,axis=2)
v_monthlymean =np.nanmean(v_mean,axis=2)

##########读取经纬度############################################################


##########计算月均的实际风速####################################################
w =  np.sqrt(u_monthlymean*u_monthlymean+v_monthlymean*v_monthlymean)
##########封装绘图函数、保护绘图区域、投影、陆地、标题、字体大小#########################
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='white'
                                        ,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]+10, xstep),crs=ccrs.PlateCarree())
    ax.set_yticks(np.arange(box[2], box[3]+10, 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=12,
                    bottom=True, left=True, right=True, top=True)
    return ax

##########绘图区域、投影陆地的分辨率、经度、纬度的间隔############################
box = [130, 180, -20, 20]  
scale = '50m'            
xstep, ystep = 10, 10 




##########进行绘图，设置一个画板##################################

fig =plt.figure(figsize=(12,10),dpi=100)

ax=fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=180))

make_map(ax,'monthly_mean-surface winds')

ax.set_xlabel('Longitude',fontsize=14)  # 添加经度标签
ax.set_ylabel('Latitude',fontsize=14)   # 添加纬度标签


step=10# 设置风场的数据读取间隔

c=ax.contourf(lon,lat,w,transform=ccrs.PlateCarree(),extend='both')  #绘制填色图
    
ax.quiver(x[::step,::step],y[::step,::step],u_monthlymean[::step,::step],v_monthlymean[::step,::step],pivot='mid',transform=ccrs.PlateCarree(),
          width=0.002,scale=200,headlength=5,headwidth=4,
        )  #绘制风场箭头
    
#添加colorbar
cb=fig.colorbar(c,shrink=0.75,pad=0.04)   
    
cb.ax.tick_params(labelsize=10)  #设置colorbar的字体大小
cb.ax.set_title('$m/s$',fontsize=15)  #设置colorbar的标题，字体大小