一、前言
第五代再分析数据集(ERA5)提供4类土壤湿度变量,变量名为“Volumetric soil water layer”
包括4个深度层,分别为0~7、7~28、28~100、100~289 cm,单位为土壤体积含水量(m3/m3)
那么如果我需要0~100 cm的土壤湿度数据,应该怎么办呢?
根据文献阅读,方法其实非常简单,只需将土层厚度作为权重因子再求和即可
数据下载链接如下,我这里计算0~100 cm的土壤湿度
二、计算方法
三、Python脚本
import os
import rasterio
import numpy as np
# 定义输入文件夹路径
folder1 = r'输入路径' #Volume of water in soil layer 1 (0-7 cm)
folder2 = r'输入路径' #Volume of water in soil layer 2 (7-28 cm)
folder3 = r'输入路径' #Volume of water in soil layer 3 (28-100 cm)
# 定义输出文件夹路径
output_folder = r'输出路径'
# 获取文件列表
files1 = sorted([f for f in os.listdir(folder1) if f.endswith('.tif')])
files2 = sorted([f for f in os.listdir(folder2) if f.endswith('.tif')])
files3 = sorted([f for f in os.listdir(folder3) if f.endswith('.tif')])
# 确保每个文件夹中的文件数量相同
assert len(files1) == len(files2) == len(files3), "文件夹中的文件数量不一致"
# 批量处理每组栅格文件
for file1, file2, file3 in zip(files1, files2, files3):
path1 = os.path.join(folder1, file1)
path2 = os.path.join(folder2, file2)
path3 = os.path.join(folder3, file3)
# 打开栅格文件
with rasterio.open(path1) as src1, rasterio.open(path2) as src2, rasterio.open(path3) as src3:
assert src1.shape == src2.shape == src3.shape, "栅格文件的形状不一致"
# 读取栅格数据
data1 = src1.read(1)
data2 = src2.read(1)
data3 = src3.read(1)
# 逐格点加权相加
sum_data = (7*data1 + 21*data2 + 72*data3) / 100
# 更新栅格元数据
profile = src1.profile
profile.update(dtype=rasterio.float32)
# 定义输出文件路径
output_path = os.path.join(output_folder, file1)
# 写入新的栅格文件
with rasterio.open(output_path, 'w', **profile) as dst:
dst.write(sum_data.astype(rasterio.float32), 1)
print("栅格运算完成,文件已输出到:", output_folder)
结果:
四、参考
[1] Z.Q. Zhong, B. He, Y.P. Wang, et al, Disentangling the effects of vapor pressure deficit on northern terrestrial vegetation productivity.Sci. Adv.9,eadf3166(2023).
[2] 金令, 王永芳, 郭恩亮等, ERA5、GLDAS和FLDAS土壤湿度数据在内蒙古地区的适用性评价. 节水灌溉, 2023(3): 53-60.