【SEBAL模型】地面特征参数反演01--大气透射率

1. 概念与原理

参考文献：
阮俊. 基于SEBAL模型的西北内陆地区蒸散量研究[D].华中科技大学,2013.

大气透射率是指透过大气层后的辐射强度和入射前的辐射强度之比。在通常情况下大气单向透射率约为0.55~0.85，对于海拔较低且面积较小时可以取值0.75，对海拔起伏很大，则可以通过数字高程计算大气透射率。

2. 数据准备

只需准备研究区DEM数据即可，本实验将DEM作为观测点海拔高度

3. 代码实现

因为只有一个公式，实现起来其实也较为简单。

#tau 𝜏_𝑆𝑊=0.75+2×10^(−5) 𝑧
#𝑧 ：测量点的海拔高度
import os
import numpy as np
from osgeo import gdal
import glob
import time
import math
import copy

#1. 读取DEM值
path = r"D:\DATA\SEBAL\InitialData\DEM.tif"  #DEM影像
dem=gdal.Open(path)  #打开

#读取行列数和投影
col = dem.RasterXSize
row = dem.RasterYSize
geotrans = dem.GetGeoTransform()
proj = dem.GetProjection() #地图投影

#将gdal.dataset转换为array 即将tif影像转换为数组
dem_value= np.array(dem.ReadAsArray(0, 0, col, row), dtype=float)  # 将数据写成数组，对应栅格矩阵

#2. 复制文件的DEM值
tau=copy.deepcopy(dem_value)  #只复制数值，不复制地址
x=len(tau)   #一维长度
y=len(tau[0])  #二维长度

#3. 计算
for i in range(x):
    for j in range(y):
        tau[i][j]=0.75+2*math.pow(10, -5)*dem_value[i][j]

#4. 写入新的影像文件
filename_out=r"D:\DATA\SEBAL\DaQiTouSheLv\tau.tif"
driver=gdal.GetDriverByName('GTiff') #将GTiff图纸实体化为driver，括号内为文件类型
 
dataset_out=driver.Create(filename_out,col,row,1,gdal.GDT_Float32) #这里的1指的是创建的栅格数据文件只能容纳一个波段
#定义仿射矩阵和投影信息
dataset_out.SetGeoTransform(geotrans)
dataset_out.SetProjection(proj)
#写入像素值
band_out=dataset_out.GetRasterBand(1)
band_out.WriteArray(tau)  #把数组写入到栅格数据相应的波段中
#清空缓存，关闭数据集
dem.FlushCache()
dataset_out.FlushCache() #把内存中的数据写入到文件中

print("已成功！")