写在前面
在进行一些评估工作时，需要用到可降水量这个数据，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]}')
ax.set_aspect(2)
cbar = fig.colorbar(contour, ax=axes, orientation='horizontal', pad=0.1, shrink=0.9, aspect=60)
cbar.set_label('kg kg**-1')
# plt.show()
感觉和自己算的还是有的类似的，有意思的是对比ERA5数据，可以明显的发现，cmip6这个三个model在陆地上的模拟都存在明显的差别，难以模拟出陆地上的可降水量，在海上赤道中东太平洋的prw信号模拟也较弱。
与ERA5-TCWV的空间分布以及纬向平均结果是比较一致，可以用来与cmip6对比
读取下载的TCWV数据，并选择同样的时间，插值相同的分辨率与计算得到的结果进行对比
era5_path = r'I://download.nc'
era5_tcwv = xr.open_dataset(era5_path).sortby('latitude').tcwv.sel(time='2010',
latitude=slice(-40, 42))
era5_tcwv_lon = convert_lon(era5_tcwv).interp(longitude=np.arange(0, 360, 2),latitude=np.arange(-40, 42, 2))
era5_tcwv_lon
print(era5_tcwv_lon.min(), era5_tcwv_lon.max())
dataset = [era5_tcwv_lon[0], sph_trap[0], era5_tcwv_lon[0] - sph_trap[0]]
data_title = ['ERA-5-TCWV', 'ERA-5-SPH', 'diff']
fig, axes = plt.subplots(1, 3, subplot_kw={'projection': ccrs.PlateCarree(180)}, dpi=200, figsize=(15, 6))
# 绘制前两个子图并共用一个颜色条
for i, (data, ax) in enumerate(zip(dataset[:2], axes[:2])):
contour = data.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'{data_title[i]}')
ax.set_aspect(2)
cbar = fig.colorbar(contour, ax=axes[:2], orientation='horizontal', pad=0.1, shrink=0.9, aspect=60)
cbar.set_label('kg kg**-1')
contour_diff = dataset[2].plot.contourf(ax=axes[2],
levels=np.linspace(-4, 4,21),
cmap='RdBu_r', extend='both', add_colorbar=False, transform=ccrs.PlateCarree())
axes[2].coastlines()
gl = axes[2].gridlines(draw_labels=True, linestyle='--') # 添加经纬度网格并显示标签
gl.right_labels = False
gl.top_labels = False
axes[2].set_title(f'{data_title[2]}')
axes[2].set_aspect(2)
cbar_diff = fig.colorbar(contour_diff, ax=axes[2], orientation='horizontal', pad=0.1, shrink=0.9, aspect=40)
cbar_diff.set_label('kg kg**-1')
plt.show()
可以发现的是，通过比湿计算得到的可降水量与TCWV在陆地上主要的偏差也在陆地上，很奇怪。在还要上的偏差倒是比较小
对于纬向平均和经向评价结果，也能得到类似的模拟情况，在纬向平均的35°E、280°E附近存在明显的差异，分布对应于南美和南非大陆，很有意思
dataset = [era5_tcwv_lon[0], sph_trap[0], era5_tcwv_lon[0] - sph_trap[0]]
data_title = ['ERA-5-TCWV', 'ERA-5-SPH', '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 = data.mean('latitude').plot(ax=ax, label=data_title[i], color=colors[i])
ax.set_title('纬向平均值')
ax.set_xlabel('经度')
ax.set_ylabel('值')
ax.legend()
plt.show()
总结
总的来说，通过从xarray和numpy的角度，基于比湿数据得到的可降水量数据与ERA5中的TCWV的空间分布以及纬向平均结果是比较一致，但是在陆地上存在明显的差异。当然， 这里测试的数据只是2010年的1月的数据，没有从更长时间范围上进行评估，只是简单验证一下计算方法的可靠性。 在后续评估中，感觉可以用TCWV来与cmip6的Prw进行背景场的评估。
相关数据和测试代码放到GitHub上
https://github.com/Blissful-Jasper/jianpu_record
★https://confluence.ecmwf.int/pages/viewpage.action?pageId=56667102
给我点在看的人
越来越好看
