《地表温度反演idl程序》由会员分享，可在线阅读，更多相关《地表温度反演idl程序（5页珍藏版）》请在金锄头文库上搜索。
1、PRO LST ; 从文件夹中读取输入数据的文件名，数据顺序为 2、4、6. imagenames=findfile(E:rsld2热岛*, count=count) outpath = E:rsld2200802sz_lstfor k =3,count-1,2 do begin print,imagenamesk; 得到图像的行列数目及投影信息 envi_open_file,imagenamesk,R_FID=image_fid, /no_realize if (image_fid eq -1) then return envi_file_query, image_fid, ns=ns, n
2、l=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
3、_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.
4、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 = flta
5、rr(ns,nl) t10_32 = fltarr(ns,nl)if index_w10 ne -1 then begin t10_31index_w1 = 0.99513-0.08082*windex_w1 t10_32index_w1 = 0.99376-0.11369*windex_w1 endifif index_w20 ne -1 then begin t10_31index_w2 = 1.08692-0.12759*windex_w2 t10_32index_w2 = 1.07900-0.15925*windex_w2 endifif index_w30 ne -1 then be
6、gin t10_31index_w3 = 1.07268-0.12571*windex_w3 t10_32index_w3 = 0.93821-0.12613*windex_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 = fltar
7、r(ns,nl)if index_31_t10 ne -1 then begin dt_31index_31_t1 = 0.08 endifif index_31_t20 ne -1 then begin ;下面的公式中的 0.000325 原公式中为 0.00325，因本计算数据中亮温为实际值的 10 倍，所以多乘一个 0.1 dt_31index_31_t2 = -0.05+0.000325*(data_31index_31_t2-2780) endifif index_31_t30 ne -1 then begin dt_31index_31_t3 = -0.05 endifindex_
8、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_t10 ne -1 then begin dt_32index_32_t1 = 0.095 endifif index_32_t20 ne -1 then begin ;下面的公式中的 0.0004 原公式中为 0.004，因本计算数据中亮温为实际值的10 倍，所以多乘一个 0.1
9、 dt_32index_32_t2 = -0.065+0.0004*(data_32index_32_t2-2780) endifif index_32_t30 ne -1 then begin dt_32index_32_t3 = -0.065 endif;大气透过率的传感器视角校正函数 angle = 20 dt = -0.00322+(3.0967*(10.0(-5)*(angle2);大气透过率 t_31 = t10_31+dt_31-dt t_32 = t10_32+dt_32-dt;计算 ndvi ndvi = fltarr(ns,nl) ndvi = (data_2-data_1
10、)*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*(
11、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_31index_water = 0.99683 e_32index_water = 0.99254 endif;计算地表温度 C_31 = e_31*t_31 C_32 = e_32*t_32D_31 = (1-t_31)*(1+(1-e_31)*t_31) D_32 = (1-t_32)*(1+
12、(1-e_32)*t_32)a_31 = -64.60363 b_31 = 0.440817 a_32 = -68.72575 b_32 = 0.473453A0 = 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_3
13、1*(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_bad10 ne -1 then begin Tindex_bad1 = 99 endifif index_bad20 ne -1 then begin Tindex_bad2 = 99 endifif index_water ne -1 then begin Tindex_water = 99 endifout_name=E:rsld220080422lstENVI_WRITE_ENVI_FILE,t,out_name=out_name,map_info=map_info,R_FID = f envi_file_mng, id=f, /removeendforprint,-end- end
