地表温度反演IDL程序

PRO LST

; 从文件夹中读取输入数据的文件名，数据顺序为2、4、6..... imagenames=findfile('E:\rs\ld2\热岛\*', count=count)

outpath = 'E:\rs\ld2\200802\sz_lst\'

for k =3,count-1,2 do begin

print,imagenames[k]

; 得到图像的行列数目及投影信息

envi_open_file,imagenames[k],R_FID=image_fid, /no_realize

if (image_fid eq -1) then return

envi_file_query, image_fid, ns=ns, nl=nl

map_info = envi_get_map_info(fid=image_fid)

dims = [-1, 0, ns-1, 0, nl-1]

;1,2,19波段的反射率

data_1 = envi_get_data(fid=image_fid,dims=dims,pos=0)

data_2 = envi_get_data(fid=image_fid,dims=dims,pos=1)

data_19 = envi_get_data(fid=image_fid,dims=dims,pos=20)

;31,32波段的亮温

data_31 = envi_get_data(fid=image_fid,dims=dims,pos=32)

data_32 = envi_get_data(fid=image_fid,dims=dims,pos=33)

index_bad1= where((data_1 eq 65534) or (data_2 eq 65534) or (data_19 eq 65534) or $

(data_31 eq 65534) or (data_32 eq 65534))

envi_file_mng, id=image_fid, /remove

;计算大气水含量

;w=((alfa-ln(ref19/ref2))/beta)2 2次幂

;alfa=0.02 beta=0.651

a=data_19*1.0/data_2

b=ALOG(a)

w=((0.02-b)/0.651)^2

;计算大气透过率

;传感器视角为10度的星下大气透过率，在水汽含量为0.4-2.0，2.0-4.0，4.0-6.0的计算方程

index_w1 = where((w ge 0.4) and (w lt 2.0))

index_w2 = where((w ge 2.0) and (w lt 4.0))

index_w3 = where((w ge 4.0) and (w le 6.0))

index_bad2 = where((w lt 0.4) or (w gt 6.0))

t10_31 = fltarr(ns,nl)

t10_32 = fltarr(ns,nl)

if index_w1[0] ne -1 then begin

t10_31[index_w1] = 0.99513-0.08082*w[index_w1]

t10_32[index_w1] = 0.99376-0.11369*w[index_w1]

endif

if index_w2[0] ne -1 then begin

t10_31[index_w2] = 1.08692-0.12759*w[index_w2]

t10_32[index_w2] = 1.07900-0.15925*w[index_w2]

endif

if index_w3[0] ne -1 then begin

t10_31[index_w3] = 1.07268-0.12571*w[index_w3]

t10_32[index_w3] = 0.93821-0.12613*w[index_w3]

endif

;大气透过率温度校正函数

;因数据中亮温为实际亮温的10倍，阈值都乘以10倍，如318k在计算中用3180 index_31_t1 = where(data_31 gt 3180)

index_31_t2 = where((data_31 le 3180) and (data_31 ge 2780))

index_31_t3 = where(data_31 lt 2780)

dt_31 = fltarr(ns,nl)

if index_31_t1[0] ne -1 then begin

dt_31[index_31_t1] = 0.08

endif

if index_31_t2[0] ne -1 then begin

;下面的公式中的0.000325原公式中为0.00325，因本计算数据中亮温为实际值的10倍，所以多乘一个0.1

dt_31[index_31_t2] = -0.05+0.000325*(data_31[index_31_t2]-2780)

endif

if index_31_t3[0] ne -1 then begin

dt_31[index_31_t3] = -0.05

endif

index_32_t1 = where(data_32 gt 3180)

index_32_t2 = where((data_32 le 3180) and (data_31 ge 2780))

index_32_t3 = where(data_32 lt 2780)

dt_32 = fltarr(ns,nl)

if index_32_t1[0] ne -1 then begin

dt_32[index_32_t1] = 0.095

endif

if index_32_t2[0] ne -1 then begin

;下面的公式中的0.0004原公式中为0.004，因本计算数据中亮温为实际值的10倍，所以多乘一个0.1

dt_32[index_32_t2] = -0.065+0.0004*(data_32[index_32_t2]-2780)

endif

if index_32_t3[0] ne -1 then begin

dt_32[index_32_t3] = -0.065

endif

;大气透过率的传感器视角校正函数

angle = 20

dt = -0.00322+(3.0967*(10.0^(-5))*(angle^2))

;大气透过率

t_31 = t10_31+dt_31-dt

t_32 = t10_32+dt_32-dt