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
;计算ndvi ndvi = fltarr(ns,nl)
ndvi = (data_2-data_1)*1.0/(data_2+data_1) ;地表比辐射率计算要区分水陆像元 index_water = where(ndvi lt -0.01) ;计算陆地像元的比辐射率 ;计算植被覆盖率 p = (ndvi-0.15)/(0.9-0.15) ;比辐射率
;(p*(p le (1-p))+(1-p)*(p gt (1-p))) 即 min(p,1-p) e_31 = p*(0.92762+0.07033*p)*0.98672+$ (1-p)*(0.99782+0.08362*p)*0.96767+$ 0.003796*(p*(p le (1-p))+(1-p)*(p gt (1-p))) e_32 = p*(0.92762+0.07033*p)*0.98990+$ (1-p)*(0.99782+0.08362*p)*0.97790+$ 0.003796*(p*(p le (1-p))+(1-p)*(p gt (1-p))) ;加入水体的比辐射率
if index_water ne -1 then begin e_31[index_water] = 0.99683 e_32[index_water] = 0.99254 endif
;计算地表温度 C_31 = e_31*t_31 C_32 = e_32*t_32
D_31 = (1-t_31)*(1+(1-e_31)*t_31) D_32 = (1-t_32)*(1+(1-e_32)*t_32) a_31 = -64.60363 b_31 = 0.440817 a_32 = -68.72575 b_32 = 0.473453
A0 = a_31*D_32*(1-C_31-D_31)/(D_32*C_31-D_31*C_32)-$ a_32*D_31*(1-C_32-D_32)/(D_32*C_31-D_31*C_32) A1 = 1+D_31/(D_32*C_31-D_31*C_32)+$
b_31*D_32*(1-C_31-D_31)/(D_32*C_31-D_31*C_32) A2 = D_31/(D_32*C_31-D_31*C_32)+$
b_32*D_31*(1-C_32-D_32)/(D_32*C_31-D_31*C_32)
;原公式为T = A0+A1*data_31-A2*data_32，但因数据中亮温为实际亮温的10倍
;所以要除以10。最后减去273转换为摄氏温度 T = A0+A1*data_31*0.1-A2*data_32*0.1-273 ;标定坏数据为99
if index_bad1[0] ne -1 then begin T[index_bad1] = 99 endif
if index_bad2[0] ne -1 then begin T[index_bad2] = 99 endif
if index_water ne -1 then begin T[index_water] = 99 endif
out_name='E:\\rs\\ld2\\20080422lst'
ENVI_WRITE_ENVI_FILE,t,out_name=out_name,map_info=map_info,R_FID = f
envi_file_mng, id=f, /remove endfor
print,'------------end---------------' end