Python | 计算可降水量

写在前面

在进行一些评估工作时，需要用到可降水量这个数据，Cmip6的模式中是可以输出的，即 prw 这个变量，定义为

Water Vapor Path [kg m-2]

但是，貌似ERA5-monthly数据中是没有直接叫可降水量的数据，以下列出了相关单位为[kg ]的变量

搜了一下相关资料，感觉只有 TCWV (total column of water vapour), 整层气柱水汽总量； TCW (total column of water), 整层气柱水含量； VIWV (vertical integral of water vapour), 水汽垂直积分比较相关，其中，变量 TCWV (水汽总量)和 VIWV (水汽垂直积分)基本相同，但由不同的软件生成。

在ERA5-monthly-single level上是可以下载total column of water vapour数据以及比湿( specific_humidity )数据的，而比湿是可以计算得到可降水量的；所以，这里通过手动从比湿计算得到可降水量数据，与TCWV数据和cmip6中的一些模式数据进行对比，来验证手动计算可降水量的结果。

比湿计算可降水量公式

可降水量，定义为垂直积分的比湿(q)，主要通过

• 转换气压层单位为pa
• 使用梯形积分计算整层水汽积分

其中，单位换算如下：

• q 是比湿（单位为 kg kg
• 是地表气压
• 是顶层气压
• g为重力加速度，单位：m/s

其中，pa单位可以换算为

= =

读取数据

先简单读取比湿数据，并转换经度从-180~180到0~360，将纬度从南到北排序

import xarray as xr
import numpy  as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

sph = xr.open_dataset("I://sph.monthly.2010_2019.nc").q
def convert_lon(ds,lons='longitude',lats='latitude'):
    ds = ds.sortby(lats)
    lon_name = lons  #你的nc文件中经度的命名
    ds['longitude_adjusted'] = xr.where(ds[lon_name] < 0, ds[lon_name]%360,\
                        ds[lon_name])
    ds = (
    ds
    .swap_dims({lon_name: 'longitude_adjusted'})
    .sel(**{'longitude_adjusted': sorted(ds.longitude_adjusted)})
    .drop_vars(lon_name))
    ds = ds.rename({'longitude_adjusted': lon_name})

    return ds

sph_convert = convert_lon(sph)
sph_convert

• 空间插值

为了方便后续计算效率，将数据插值到2x2分辨率，并看一下插值的空间分布图

sph_interp = sph_convert.sel(level=slice(100,1000)).interp(longitude=np.arange(0,360,2),latitude=np.arange(-40,42,2))
sph_interp 
fig, ax = plt.subplots(subplot_kw={'projection': ccrs.PlateCarree(180)}, dpi=200,figsize=(8, 6))
contour = sph_interp.sel(level=850)[0].plot.contourf(ax=ax, levels=11, cmap='RdBu_r', extend='both', add_colorbar=False,
                                              transform=ccrs.PlateCarree())
ax.coastlines()
gl = ax.gridlines(draw_labels=True




    
, linestyle='--')  # 添加经纬度网格并显示标签
gl.right_labels = False  # 关闭右侧标签
gl.top_labels = False    # 关闭顶部标签

# 添加横轴方向的colorbar
cbar = fig.colorbar(contour, ax=ax, orientation='horizontal', pad=0.1, shrink=0.9,aspect=40)
cbar.set_label('kg kg**-1')  
plt.show()

下面提供两种计算方法， 一种是xarray的，一种是numpy的，双重验证

xarray

def xarray_intergrated(sph):

    pressure = (sph.level * 100) # 转换为 Pa
    print(pressure)
    # 定义重力加速度常数
    g = 9.81  # m/s²

    # 计算dp（差分），并取绝对值
    dp = pressure.diff(dim='level')
   
    # 计算垂直积分 (Column-integrated specific humidity)
    column_integrated_sph = (sph * dp).sum(dim='level') / g

    # 输出结果，单位为 kg/m²
    print(column_integrated_sph.min(), column_integrated_sph.max())
    return column_integrated_sph
sph_trap = xarray_intergrated(sph_interp)
sph_trap

简单绘图看看一下空间分布

fig, ax = plt.subplots(subplot_kw={'projection': ccrs.PlateCarree(180)}, dpi=200,figsize=(8, 6))

contour = sph_trap[0].plot.contourf(ax=ax, levels=np.linspace(20,70,51), cmap='RdBu_r', extend='both', add_colorbar=False,
                                              transform=ccrs.PlateCarree())
ax.coastlines()
gl = ax.gridlines(draw_labels=True, linestyle='--')  # 添加经纬度网格并显示标签
gl.right_labels = False  # 关闭右侧标签
gl.top_labels = False    # 关闭顶部标签
ax.set_aspect(2)
# 添加横轴方向的colorbar
cbar = fig.colorbar(contour, ax=ax, orientation='horizontal', pad=0.1, shrink=0.9,aspect=40)
cbar.set_label('kg kg**-1')  

plt.show()

numpy

def numpy_intergrated(sph_interp):
    
    pressure_levels = sph_interp.level * 100  # 将 millibars 转换为 Pa

    # 重力加速度
    g = 9.81  # m/s²

    # 使用梯形积分法计算整层水汽积分
    sph_trap = np.trapz(sph_interp, pressure_levels, axis=1) / g
    print(sph_trap.min(),sph_trap.max())
    return sph_trap

sph_trap_numpy = numpy_intergrated(sph_interp)    

同样看一下空间分布

fig, ax = plt.subplots(subplot_kw={'projection': ccrs.PlateCarree(180)}, dpi=200,figsize=(8, 6))

contour = ax.contourf(sph_trap_numpy[0], 
                      levels=np.linspace(20,70,51), 
                      
                      cmap='RdBu_r', extend='both', 
                                              transform=ccrs.PlateCarree())
ax.coastlines()
gl = ax.gridlines(draw_labels=True, linestyle='--')  # 添加经纬度网格并显示标签
gl.right_labels = False  # 关闭右侧标签
gl.top_labels = False    # 关闭顶部标签
ax.set_aspect(2)
# 添加横轴方向的colorbar
cbar = fig.colorbar(contour, ax=ax, orientation='horizontal', pad=0.1, shrink=0.9,aspect=40)
cbar.set_label('kg kg**-1')  

plt.show()

从纬向平均的情况来看看是否有明显的差异

dataset = [sph_trap[0].mean('latitude'), np.nanmean(sph_trap_numpy[0],axis=0), sph_trap[0].mean('latitude')- np.nanmean(sph_trap_numpy[0],axis=0)]
data_title = ['Prw-Xarray', 'Prw-Numpy', 'diff']

colors = ['b', 'g', 'r']  

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
fig, ax = plt.subplots(dpi=200, figsize=(10, 6))

# 绘制三个数据的纬向平均值在同一个图上
for i, data in enumerate(dataset):
    line = ax.plot(data, label=data_title[i], color=colors[i])

ax.set_title('纬向平均值')
ax.set_xlabel('经度')
ax.set_ylabel('值')
ax.legend()

plt.show()

可以发现，总体上趋势还是一致的，基于Xarray的计算方法稍微比基于Numpy的结果要大一点，不影响总体结论

与CMIP6的空间patter对比

下面，找了几个CMIP6的prw月平均资料，选择同样的时间，同样的分辨率来看看空间分布

import xarray as xr
import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs

model_paths = [
    "I://EC-Earth3_historical_r1i1p1f1_1997-2014_interpolated.nc",
    "I://SAM0-UNICON_historical_r1i1p1f1_1997-2014_interpolated.nc",
    "I://UKESM1-0-LL_historical_r1i1p1f2_1997-2014_interpolated.nc",
]

model_names = ['EC-Earth3', 'SAM0-UNICON', 'UKESM1-0-LL']
# 读取所有模型数据
models = [xr.open_dataset(path).prw.sel(time='2010', lat=slice(-40, 40), lon=slice(0, 360))[0] for path in model_paths]
print(models[0].min(),models[0].max())

fig, axes = plt.subplots(1, 3, subplot_kw={'projection': ccrs.PlateCarree(180)}, dpi=200, figsize=(15, 6))

for i, (model, ax) in enumerate(zip(models, axes)):

    model.plot.contourf(ax=ax, levels=np.linspace(20, 70, 51), cmap='RdBu_r', extend='both',add_colorbar=False,transform=ccrs.PlateCarree())
    ax.coastlines()
    gl = ax.gridlines(draw_labels=True, linestyle='--')  # 添加经纬度网格并显示标签
    gl.right_labels = False  # 关闭右侧标签
    gl.top_labels = False    # 关闭顶部标签
    ax.set_title(f'{model_names[i]}