scipy 如何在Python中计算水汽通量收敛的垂直积分?

qvk1mo1f  于 2022-11-10  发布在  Python
关注(0)|答案(1)|浏览(607)

我是python的新手,我使用metpy.calc函数HMC_LE = (np.array(mpcalc.divergence(uq, vq, dx=dx, dy=dy)))计算了印度850 hPa气压水平下的水汽通量散度[MFD](即{d(qu)/dx}+{d(qv)/dy},q是比湿度,u、v是纬向和经向风分量)。我遵循了以下步骤:How to calculate moisture flux divergence in python。我的代码-

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cartopy.util import add_cyclic_point
import pandas as pd
import cartopy.feature as cf
from netCDF4 import Dataset 
import metpy.calc as mpcalc

ds1 = xr.open_dataset('Downloads/uwnd.mon.mean.nc')
ds2 = xr.open_dataset('Downloads/vwnd.mon.mean.nc')
ds3 = xr.open_dataset('Downloads/shum.mon.mean.nc')
lons = ds1['lon'][:]  
lats = ds1['lat'][:] 
levs = ds1['level'][2]

# --- select all JJAS months

uJJAS_ = ds1['uwnd'].sel(time=np.in1d(ds1['time.month'], [6, 7, 8,9]))
uJJAS_
uJJAS=uJJAS_[:,2,:,:]
uJJAS

# Select JJAS 850 hPa from 1981-2010

u_1=uJJAS[132:252,:,:]
u_1

# Prepare the climatology

u_1_mean= np.mean(u_1,axis=0)
u_1_mean

# --- For v

vJJAS_ = ds2['vwnd'].sel(time=np.in1d(ds2['time.month'], [6, 7, 8,9]))
vJJAS_
vJJAS=vJJAS_[:,2,:,:]
vJJAS    
v_1=vJJAS[132:252,:,:]
v_1
v_1_mean= np.mean(v_1,axis=0)
v_1_mean

# --- For specific humidity

shJJAS_ = ds3['shum'].sel(time=np.in1d(ds3['time.month'], [6, 7, 8,9]))
shJJAS_
shJJAS=shJJAS_[:,2,:,:]
shJJAS    
sh_1=shJJAS[132:252,:,:]
sh_1
sh_1_mean= np.mean(sh_1,axis=0)
sh_1_mean

uq=(u_1_mean)*(sh_1_mean)
vq=(v_1_mean)*(sh_1_mean)

# Compute dx and dy spacing for use in divergence calculation

dx, dy = mpcalc.lat_lon_grid_deltas(lons, lats)

# Calculate moisture flux divergence

HMC_LE = (np.array(mpcalc.divergence(uq, vq, dx=dx, dy=dy)))

lat=lats.to_numpy()
lon=lons.to_numpy()
u_=uq.to_numpy()
v_=vq.to_numpy()
ax2 = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
clevs = [-0.06,-0.05,-0.04,-0.03,-0.02,-0.01, 0,0.01,0.02,0.03,0.04,0.05,0.06]
vp_fill = ax2.contourf(lons,lats,HMC_LE*1000,clevs,transform=ccrs.PlateCarree(),cmap=plt.cm.RdBu_r,
extend='both')
plt.colorbar(vp_fill, orientation='vertical')
q2=ax2.quiver(lon, lat, u_, v_,width=0.003,headwidth=5, scale_units='xy',scale=25,transform=ccrs.PlateCarree())
ax2.coastlines(alpha=0.8)
ax2.set_xticks([60,70,80,90,100], crs=ccrs.PlateCarree())
ax2.set_yticks([-10,0,10,20,30,40], crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True,number_format='.0f')
lat_formatter = LatitudeFormatter()
ax2.xaxis.set_major_formatter(lon_formatter)
ax2.yaxis.set_major_formatter(lat_formatter)
ax2.set_extent([57,101,-10,45])
ax2.add_feature(cf.BORDERS,edgecolor='black',alpha=0.8)
plt.title('Moisture flux divergence at 850 hPa ', fontsize=16)
plt.show()

所得到的曲线如图所示-x1c 0d1x I具有1000百帕、925百帕、850百帕、700百帕、600百帕、500百帕、400百帕的值,300 hPa。现在我想对MFD值进行垂直积分[从1000 hPa到300 hPa气压水平]。请告诉我如何进行积分?(({d(qu)/dx}+{d(qv)/dy})dP从300百帕到1000百帕的积分。)
谢谢你,谢谢你

7z5jn7bk

7z5jn7bk1#

如果将所有级别组合到一个数组中,则应该能够使用numpy.trapzscipy.integrate中的一些更复杂的积分方法进行积分。

相关问题