区域净水汽收支计算并绘图

最新推荐文章于 2024-05-22 17:37:52 发布

豆子_Yan

最新推荐文章于 2024-05-22 17:37:52 发布

阅读量1.5k

点赞数 27

文章标签： python matplotlib

本文链接：https://blog.csdn.net/weixin_64067690/article/details/135746521

版权

文章介绍了如何使用Python和MetPy库处理ERA5数据，计算径向和纬向的净水汽通量，然后通过Cartopy进行区域划分并绘制出净水汽通量的平面分布图，同时展示了风场矢量图。

摘要由CSDN通过智能技术生成

1.净水汽通量计算理论

径向整层水汽通量： $Q_{v}=\frac{1}{g}\int_{p_{s}}^{p_{t}}qvdp$

纬向整层水汽通量： $Q_{u}=\frac{1}{g}\int_{p_{s}}^{p_{t}}qdp$

纬向西边界净水汽通量： $F_{w}=\int_{\varphi _s}^{\varphi _n}Q_{uW}Rd\varphi$

纬向东边界净水汽通量： $F_{e}=\int_{\varphi _s}^{\varphi _n}Q_{uE}Rd\varphi$

经向南边界净水汽通量： $F_{s}=\int_{\lambda _w}^{\lambda _e}Q_{vS}R\cos \varphi _{s}d\lambda$

经向北边界净水汽通量： $F_{n}=\int_{\lambda _w}^{\lambda _e}Q_{vN}R\cos \varphi _{n}d\lambda$

其中，为地球半径， $\varphi$ 为纬度， $\lambda$ 为经度

区域净水汽通量: $F=F_{w}-F_{e}+F_{s}-F_{n}$

2.python实现以及区域分布绘图

import cartopy.crs as ccrs
import cartopy.feature as cfeature

import metpy.constants as constants
from matplotlib.ticker import MultipleLocator, FormatStrFormatter
from mpl_toolkits.basemap import Basemap
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

import numpy as np
from matplotlib import pyplot as plt
import pandas as pd
import datetime as dt
import xarray as xr
import metpy.calc as mpcalc
from metpy.units import units
plt.rcParams['font.family'] = 'Times New Roman'
fontdict = {"size":8,"weight":"bold"}

ds=xr.open_dataset(r'*/era5_u_v_w_t_rh_q2.nc') #era5
area_f=[] #区域上边界
lev = ds.level[:] # 读取气压层，单位为mb，即hPa，一维的14.
lat_b=22
lat_t=42
lon_left=70
lon_right=110
time_left="2021-07-25T00:00:00"
time_right="2021-07-26T00:00:00"
time1=time_left[8:10]
time2=time_right[8:10]
print(time1,time2)
#区域经纬度
for lat in np.arange(lat_b,lat_t+0.1,0.25):
for lon in np.arange(lon_left,lon_right+0.1,0.25):
u = ds.u.loc[time_left:time_right,:,lat+0.25:lat,lon:lon+0.25].mean(axis=0) # U风分量，单位为m/s
v = ds.v.loc[time_left:time_right,:,lat+0.25:lat,lon:lon+0.25].mean(axis=0) # V风分量，单位为m/s
q = ds.q.loc[time_left:time_right,:,lat+0.25:lat,lon:lon+0.25].mean(axis=0) # 读取比湿，单位为kg/kg
Lon=q.longitude
Lat=q.latitude
# # # 计算单层水汽通量和水汽通量散度
qv_u = u*q/(constants.g*10**-2) # g的单位为m/s2，换算为N/kg，再换算为10-2hPa·m2/kg,最终单层水汽通量的单位是kg/m•hPa•s
qv_v = v*q/(constants.g*10**-2)
total_q_u = np.trapz(qv_u,lev,axis=0) #将单位kg/(m*s)
total_q_v = np.trapz(qv_v,lev,axis=0)

F_u=-np.trapz(total_q_u,Lat,axis=0)*(6.37*10**6) #纬向积分纬度
F_v=np.trapz(total_q_v,Lon,axis=1)*(6.37*10**6)*np.cos(lat*np.pi/180)#纬向积分经度
we_f=F_u.data[0]-F_u.data[-1]+F_v.data[-1]-F_v.data[0]
area_f.append(we_f)

fontdict1 = {"size":6,"weight":"bold"}
#平均时间段内，所有高度层风场
u1 = ds.u.loc[time_left:time_right,:,lat_t:lat_b,lon_left:lon_right].mean(axis=0).mean(axis=0) # U风分量，单位为m/s
v1 = ds.v.loc[time_left:time_right,:,lat_t:lat_b,lon_left:lon_right].mean(axis=0).mean(axis=0)
ws=(u1**2+v1**2)**0.5
#水汽通量转换为2维与经纬度长度相同（纬度，经度）
a=np.array(area_f).reshape(int((lat_t-lat_b)/0.25)+1,int((lon_right-lon_left)/0.25)+1)

scale='50m'
proj = ccrs.PlateCarree()
fig = plt.figure(figsize=(5,3),dpi=600,facecolor="w")
ax1=plt.subplot(1,1,1,projection=proj)
#绘制省界
m = Basemap(llcrnrlon=lon_left,llcrnrlat=lat_b,urcrnrlon=lon_right,urcrnrlat=lat_t)
m.readshapefile(r'G:\shpfile\country1_4l_4p\Data_ipynb\bou2_4l','province',color='gray',linewidth=0.5,default_encoding='latin-1')
ax1.add_feature(cfeature.OCEAN.with_scale(scale)) #绘制海洋
ax1.add_feature(cfeature.LAND.with_scale(scale),lw=0.2,color='ghostwhite') #绘制陆地
ax1.add_feature(cfeature.COASTLINE,lw=0.25) #绘制海岸线

#绘制净水汽通量平面分布
nwf= ax1.contourf(np.arange(lon_left,lon_right+0.1,0.25),np.arange(lat_b,lat_t+0.1,0.25),a/10**6
,cmap='bwr', extend='both',levels=range(-15,16,5) )
cbar = fig.colorbar(nwf, shrink=0.7, pad=0.02,ticks=range(-15,16,5)) #colorbar设置
cbar.ax.tick_params(direction='in',length=2,labelsize=5)
cbar.set_label('Net Water Vapor Flux ( 10$^6$kg/s )', fontsize=8)
#每间隔5个数据绘制风场，目的稀疏风矢量；风速标准后使得矢量箭头大小一致
ax1.quiver(u1.longitude[::5],u1.latitude[::5], (u1/ws)[::5,::5], (v1/ws)[::5,::5], pivot='tail',width=0.002,
scale=50, color='k', headwidth=4,alpha=1,transform=ccrs.PlateCarree(),zorder=10)

ax1.set_yticks(range(lat_b,lat_t+1,5)) # 指定要显示的经纬度
ax1.set_xticks(range(lon_left,lon_right+1,5))
ax1.xaxis.set_major_formatter(LongitudeFormatter()) # 刻度格式转换为经纬度样式
ax1.yaxis.set_major_formatter(LatitudeFormatter())
ax1.tick_params(axis='both',which="major",length=2,width=1.0,color="k",direction='in',labelsize=6)
ax1.tick_params(axis='both',which="minor",length=1,width=1.0,color="k",direction='in',labelsize=6)
ax1.tick_params(top='in', right='in', which='both')
ax1.set_xlabel('Longitude',fontdict=fontdict)
ax1.set_ylabel('Latitude',fontdict=fontdict)
ax1.set_title('07-{} 00:00 UTC ~ 07-{} 00:00 UTC'.format(time1,time2),fontsize=8,loc="left" )
plt.show()

3.输出结果