Python实现水云模型批量处理

在基于微波遥感的土壤水分反演工作中，我们需要消除植被覆盖对土壤水分反演的影响，水云模型是进行此类工作的经典模型，要通过python构建水云模型消除植被覆盖影响，需要经历以下两个步骤：

一、区分植被覆盖土壤和裸露土壤

植被含水量(VWC, Vegetation Water Content)可以作为区分两类地表的重要依据，
当VWC>0时，视为植被覆盖地表，当VWC≤0时，则视为裸露地表。
𝑽𝑾𝑪=𝟏.𝟗𝟏𝟑𝟒(𝑵𝑫𝑽𝑰)^𝟐−𝟎.𝟑𝟐𝟏𝟓𝑵𝑫𝑽𝑰     (𝟎.𝟏𝟕<𝑵𝑫𝑽𝑰≤𝟎.𝟓)
𝑽𝑾𝑪=𝟒.𝟐𝟖𝟓𝟕𝑵𝑫𝑽𝑰−𝟏.𝟓𝟒𝟐𝟗     (𝑵𝑫𝑽𝑰>𝟎.𝟓)
𝑽𝑾𝑪=𝟎     (𝑵𝑫𝑽𝑰≤𝟎.𝟏𝟕)

二、根据VWC的值，对植被覆盖土壤进行水云模型的处理

以下是水云模型批量处理的实现代码：

import numpy as np
import pandas as pd
import os


def openfile(file_path):

    # 打开文件
    with open(file_path, 'r') as f:
        # 读取文件内容
        content = f.readlines()

    # 初始化空的二维数组
    two_dimensional_array = []
    # 遍历文件内容行
    for line in content[6:]:
        # 将每一行按空格分割，并将每个元素转换为浮点数
        row = [float(x) for x in line.split()]

        # 将行添加到二维数组中
        two_dimensional_array.append(row)
    # 返回二维数组
    return two_dimensional_array


def count_vwc(ndvi):
    vwc = np.zeros_like(ndvi)  # 创建一个与 ndvi 数组形状相同的数组，用于存储计算后的 VWC 值
    for i, line in enumerate(ndvi):
        for j, value in enumerate(line):
            if value <= 0.17:
                vwc[i, j] = 0
            elif 0.17 < value <= 0.5:
                vwc[i, j] = 1.9314 * pow(value, 2) - 0.3215 * value
            else:
                vwc[i, j] = 4.2857 * value - 1.5429
    return vwc


def count_tao2(vwc, theta):
    vwc_array = np.array(vwc)
    theta_array = np.array(theta)
    B = 0.091
    theta_array = np.radians(theta_array)  # 必须角度转弧度
    tao2 = np.exp(-2*B*vwc_array*(1 / np.cos(theta_array)))
    return tao2


def count_sigma_veg(vwc, tao2, theta):
    A = 0.0012
    theta_array = np.array(theta)
    theta_array = np.radians(theta_array)
    vwc_array = np.array(vwc)
    tao2_array = np.array(tao2)
    sigma_veg = A*vwc_array*np.cos(theta_array)*(1-tao2_array)
    return sigma_veg


def count_sigma_soil(vv, veg, tao2):
    vv_array = np.array(vv)
    veg_array = np.array(veg)
    tao2_array = np.array(tao2)
    sigma_soil = (vv_array - veg_array)/tao2_array
    return sigma_soil


def process_folder(input_folder_vv, input_folder_vh, input_folder_ndvi, input_folder_angle, output_folder):
    for file_name_vv in os.listdir(input_folder_vv):
        if file_name_vv.endswith('.txt'):
            file_name_vh = file_name_vv.replace('VV', 'VH')
            file_path_vv = os.path.join(input_folder_vv, file_name_vv)
            file_path_vh = os.path.join(input_folder_vh, file_name_vh)
            # 构建NDVI文件名
            file_name_ndvi = file_name_vv.replace('VV', 'NDVI')
            file_path_ndvi = os.path.join(input_folder_ndvi, file_name_ndvi)

            file_name_angle = file_name_vv.replace('VV', 'angle')
            file_path_angle = os.path.join(input_folder_angle, file_name_angle)

            VV = openfile(file_path_vv)
            VH = openfile(file_path_vh)
            ndvi = openfile(file_path_ndvi)
            angle = openfile(file_path_angle)
            VWC = count_vwc(ndvi)

            # 将 VV、VH 和 VWC 转换为一维数组
            VV = np.array(VV)
            VH = np.array(VH)
            angle = np.array(angle)
            vv_flat = VV.flatten()
            vh_flat = VH.flatten()
            vwc_flat = VWC.flatten()
            angle_flat = angle.flatten()

            # 创建一个 DataFrame 保存数据
            df = pd.DataFrame({'vv': vv_flat, 'vh': vh_flat, 'vwc': vwc_flat, 'angle':angle_flat})

            # 根据 VWC 的值判断处理方式
            if any(vwc_flat):
                TAO2 = count_tao2(VWC, angle)
                Sigma_VEG = count_sigma_veg(VWC, TAO2, angle)
                Sigma_Soil_vv = count_sigma_soil(VV, Sigma_VEG, TAO2)
                Sigma_Soil_vh = count_sigma_soil(VH, Sigma_VEG, TAO2)
                # 将计算后的数据添加到 DataFrame 中
                df['sigma_soil_vv'] = Sigma_Soil_vv.flatten()
                df['sigma_soil_vh'] = Sigma_Soil_vh.flatten()
            else:
                # 将 vwc=0 时的 vv 和 vh 直接输出
                df['sigma_soil_vv'] = vv_flat
                df['sigma_soil_vh'] = vh_flat

            # 提取输入文件的基本名称
            base_name = os.path.splitext(file_name_vv)[0]
            output_file_path = os.path.join(output_folder, f'{base_name}_output.csv')
            df.to_csv(output_file_path, index=False)


if __name__ == '__main__':
    input_folder_vv = "E:/Image/实验/ASCII/ASCII_VV2"
    input_folder_vh = "E:/Image/实验/ASCII/ASCII_VH2"
    input_folder_ndvi = "E:/Image/实验/ASCII/ASCII_ndvi2"
    input_folder_angle = "E:/Image/实验/ASCII/ASCII_angle2"
    output_folder = "E:/Image/实验/WCM3/"
    process_folder(input_folder_vv, input_folder_vh, input_folder_ndvi, input_folder_angle, output_folder)

 参考文献

Attema E P W, Ulaby F T. Vegetation modeled as a water cloud[J]. Radio Science, 1978, 13(2): 357-364.