Python——基于CRU数据的饱和水汽压差（VPD）计算及批量输出

一、定义

饱和水汽压差（Vapor Pressure Deficit，简称VPD） 是指在一定温度下，饱和水汽压与空气中的实际水汽压之间的差值。VPD对于植被生理状态具有重要影响。

简言之，VPD = 饱和水汽压 - 实际水汽压

最近在利用机器学习分析不同环境因子对碳通量、SIF的重要性，其中涉及VPD，但苦于目前没有一套长时序的全球覆盖VPD产品，所以我干脆使用CRU数据自行计算VPD。

---

二、计算方法

（一）公式

我所使用的公式是世界气象组织推荐的 Goff-Gratch 公式，需要注明的是，该公式有改良版和未改良版，其中改良版区分了temp＜0℃和temp＞0℃的情景（即水面和冰面的气压），我这里使用了未改良版的公式，因为我的研究只考虑植被覆盖地区，不考虑水体区域。

（二）数据

数据准备：① CRU temp（温度）数据      ② CRU vap（实际水汽压）数据

（上述为CRU TS v.4.03产品，该数据全球范围覆盖，各参数每月得一个平均值）

数据下载地址：CRU数据下载地址（包含温度和实际水汽压）

提醒：原数据为0.5°×0.5°/monthly分辨率，格式为nc，下载后需要转为栅格格式

---

三、Python脚本

（一）计算饱和水汽压

import os
import numpy as np
from osgeo import gdal

# Goff-Gratch公式的计算函数
def goff_gratch(t_kelvin):
    return np.power(10, 0.1079574e2 * (1-(273.16 / t_kelvin)) - 5.02800 * np.log10(t_kelvin/273.16) \
    + 0.150475e-3 * (1-np.power(10, -8.2969 * ((t_kelvin / 273.16)-1))) + 0.42873e-3 \
    * (np.power(10, 4.76955 * (1-(273.16 / t_kelvin))) - 1) + 0.78614)

def process_tif(input_path, output_path):
    # 打开输入TIFF文件
    ds = gdal.Open(input_path)
    band = ds.GetRasterBand(1)
    temp_data = band.ReadAsArray()

    # 将温度从摄氏度转换为开尔文
    temp_kelvin = temp_data + 273.15

    # 计算饱和水汽压
    es = goff_gratch(temp_kelvin)

    # 创建输出TIFF文件
    driver = gdal.GetDriverByName('GTiff')
    out_ds = driver.Create(output_path, ds.RasterXSize, ds.RasterYSize, 1, gdal.GDT_Float32)
    out_ds.SetGeoTransform(ds.GetGeoTransform())
    out_ds.SetProjection(ds.GetProjection())

    # 写入数据到输出文件
    out_band = out_ds.GetRasterBand(1)
    out_band.WriteArray(es)
    out_band.SetNoDataValue(-9999)
    out_band.FlushCache()

    # 关闭数据集
    del ds
    del out_ds

def main(input_folder, output_folder):
    if not os.path.exists(output_folder):
        os.makedirs(output_folder)

    for file_name in os.listdir(input_folder):
        if file_name.endswith('.tif'):
            input_path = os.path.join(input_folder, file_name)
            output_path = os.path.join(output_folder, f"es_{file_name}")
            process_tif(input_path, output_path)
            print(f"已处理 {file_name}")

if __name__ == "__main__":
    temp_folder = r'在此修改'  # 输入文件夹路径（温度）（可以多张tif）
    output_folder = r'在此修改'  # 输出文件夹路径（实际水汽压）（可以多张tif）
    main(input_folder, output_folder)

注意事项：Goff-Gratch公式输出为hPa，我的脚本输出的也为hPa，需要kPa的可自行修改

这里只需修改输入数据地址（即下载好的CRU温度和实际水汽压数据）

我这里计算饱和水汽压和计算VPD是两个独立的脚本，因为考虑到第二个批量作差的脚本我后续可能还会用到，所以就单独又写了一个

（二）计算VPD

import os
from osgeo import gdal

def subtract_tif(file1, file2, output_file):
    # 打开饱和水汽压和实际水汽压的文件夹
    ds1 = gdal.Open(file1)
    ds2 = gdal.Open(file2)

    band1 = ds1.GetRasterBand(1)
    band2 = ds2.GetRasterBand(1)

    data1 = band1.ReadAsArray()
    data2 = band2.ReadAsArray()

    # 计算差值（饱和 minus 实际）
    diff = data1 - data2

    # 创建输出TIFF文件
    driver = gdal.GetDriverByName('GTiff')
    out_ds = driver.Create(output_file, ds1.RasterXSize, ds1.RasterYSize, 1, gdal.GDT_Float32)
    out_ds.SetGeoTransform(ds1.GetGeoTransform())
    out_ds.SetProjection(ds1.GetProjection())

    out_band = out_ds.GetRasterBand(1)
    out_band.WriteArray(diff)
    out_band.SetNoDataValue(-9999)
    out_band.FlushCache()

    # 关闭数据集
    del ds1
    del ds2
    del out_ds

def process_folders(folder1, folder2, output_folder):
    if not os.path.exists(output_folder):
        os.makedirs(output_folder)

    files1 = sorted(os.listdir(folder1))
    files2 = sorted(os.listdir(folder2))

    for file1, file2 in zip(files1, files2):
        if file1.endswith('.tif') and file2.endswith('.tif'):
            input_path1 = os.path.join(folder1, file1)
            input_path2 = os.path.join(folder2, file2)
            output_path = os.path.join(output_folder, f"diff_{file1}")
            subtract_tif(input_path1, input_path2, output_path)
            print(f"Processed {file1} and {file2}")

if __name__ == "__main__":
    folder1 = r'在此修改'  # 饱和水汽压文件夹路径（可以多张tif）
    folder2 = r'在此修改'  # 实际水汽压文件夹路径（可以多张tif）
    output_folder = r'在此修改'  # VPD输出路径
    process_folders(folder1, folder2, output_folder)

---

参考：

1、Matlab脚本可参考另一个博主：Matlab脚本参考链接

2、饱和水汽压的计算公式参考（百度文库）：饱和水汽压计算公式介绍