地表温度反演IDL程序.docx

PRO LST; 从文件夹中读取输入数据的文件名，数据顺序为2、4、6.imagenames=findfile(E:rsld2热岛*, count=count)outpath = E:rsld2200802sz_lstfor k =3,count-1,2 do beginprint,imagenamesk; 得到图像的行列数目及投影信息envi_open_file,imagenamesk,R_FID=image_fid, /no_realizeif (image_fid eq -1) then returnenvi_file_query, image_fid, ns=ns, nl=nlmap_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.651a=data_19*1.0/data_2b=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_w10 ne -1 then begin t10_31index_w1 = 0.99513-0.08082*windex_w1 t10_32index_w1 = 0.99376-0.11369*windex_w1endifif index_w20 ne -1 then begin t10_31index_w2 = 1.08692-0.12759*windex_w2 t10_32index_w2 = 1.07900-0.15925*windex_w2endifif index_w30 ne -1 then begint10_31index_w3 = 1.07268-0.12571*windex_w3t10_32index_w3 = 0.93821-0.12613*windex_w3endif;大气透过率温度校正函数;因数据中亮温为实际亮温的10倍，阈值都乘以10倍，如318k在计算中用3180index_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_t10 ne -1 then begindt_31index_31_t1 = 0.08endifif index_31_t20 ne -1 then begin;下面的公式中的0.000325原公式中为0.00325，因本计算数据中亮温为实际值的10倍，所以多乘一个0.1dt_31index_31_t2 = -0.05+0.000325*(data_31index_31_t2-2780)endifif index_31_t30 ne -1 then begindt_31index_31_t3 = -0.05endifindex_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 begindt_32index_32_t1 = 0.095endifif index_32_t20 ne -1 then begin;下面的公式中的0.0004原公式中为0.004，因本计算数据中亮温为实际值的10倍，所以多乘一个0.1dt_32index_32_t2 = -0.065+0.0004*(data_32index_32_t2-2780)endifif index_32_t30 ne -1 then begindt_32index_32_t3 = -0.065endif;大气透过率的传感器视角校正函数angle = 20dt = -0.00322+(3.0967*(10.0(-5)*(angle2);大气透过率t_31 = t10_31+dt_31-dtt_32 = t10_32+dt_32-dt;计算ndvindvi = 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_31index_water = 0.99683 e_32index_water = 0.99254endif;计算地表温度C_31 = e_31*t_31C_32 = e_32*t_32D_31 = (1-t_31)*(1+(1-e_31)*t_31)D_32 = (1-t_32)*(1+(1-e_32)*t_32)a_31 = -64.60363b_31 = 0.440817a_32 = -68.72575b_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_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;标定坏数据为99if index_bad10 ne -1 then beginTindex_bad1 = 99endifif index_bad20 ne -1 then beginTindex_bad2 = 99end