热带森林和中纬度农田场景的完整可运行 IDL 代码

热带森林和中纬度农田场景的完整可运行IDL代码热带森林和中纬度农田场景的完整可运行IDL代码，以及ENVI批量处理适配方案，同时针对实操中高频报错（IDL读取HDF失败、ENVI批量工具闪退）提供解决方案，覆盖不同区域的全流程落地：一、热带森林MODISLST反演（完整IDL代码）场景适配：热带雨林（亚马逊/东南亚）、高水汽、高植被覆盖、多云idl;=====================================环境初始化=====================================PROTropical_Forest_LST;1.路径配置（替换为实际路径）input_path='D:/MODIS_Tropical/Forest/'output_path='D:/MODIS_Tropical/Result/'modtran_path='C:/MODTRAN6/'SETENV,'MODTRAN_DIR='+modtran_pathFILE_MKDIR,output_path;2.读取文件列表files_mod02=FILE_SEARCH(input_path,'MYD021KM*.hdf');Aqua卫星更适配热带IFN_ELEMENTS(files_mod02)EQ0THENBEGINPRINT,'未找到MYD021KM数据！'RETURNENDIF;3.日志初始化log_file=FOPEN(output_path+'Tropical_LST_Log.log',/WRITE)FPRINT,log_file,'热带森林LST反演日志|日期：'+SYSTIME()FPRINT,log_file,'==========================================';=====================================批量处理=====================================FORf_idx=0,N_ELEMENTS(files_mod02)-1DOBEGINfile_mod02=files_mod02[f_idx]file_name=FILE_BASENAME(file_mod02,'.hdf')PRINT,'处理：',file_name;-------------------------步骤1：读取HDF数据-------------------------hdf_id=HDF_OPEN(file_mod02,/READ,ERROR=err)IFerrNE0THENBEGINFPRINT,log_file,'错误：'+file_name+'-HDF文件打开失败'CONTINUEENDIF;热红外Band31/32+定标系数data31=HDFSD_GETDATA(hdf_id,'EV_1KM_Emissive',0,10)data32=HDFSD_GETDATA(hdf_id,'EV_1KM_Emissive',0,11)gain31=HDFATTR_GET(hdf_id,'EV_1KM_Emissive','gain',10)offset31=HDFATTR_GET(hdf_id,'EV_1KM_Emissive','offset',10)gain32=HDFATTR_GET(hdf_id,'EV_1KM_Emissive','gain',11)offset32=HDFATTR_GET(hdf_id,'EV_1KM_Emissive','offset',11);可见光Band1/2+定标系数data_b1=HDFSD_GETDATA(hdf_id,'EV_1KM_Reflective',0,0)data_b2=HDFSD_GETDATA(hdf_id,'EV_1KM_Reflective',0,1)gain_b1=HDFATTR_GET(hdf_id,'EV_1KM_Reflective','gain',0)offset_b1=HDFATTR_GET(hdf_id,'EV_1KM_Reflective','offset',0)gain_b2=HDFATTR_GET(hdf_id,'EV_1KM_Reflective','gain',1)offset_b2=HDFATTR_GET(hdf_id,'EV_1KM_Reflective','offset',1)HDF_CLOSE,hdf_id;-------------------------步骤2：辐射定标-------------------------L31=data31*gain31+offset31L32=data32*gain32+offset32rho_b1=data_b1*gain_b1+offset_b1rho_b2=data_b2*gain_b2+offset_b2;-------------------------步骤3：热带大气参数提取-------------------------;水汽含量（热带4~8g/cm²）file_mod07=FILE_SEARCH(input_path,'MYD07*'+STRMID(file_name,9,13)+'*.hdf')IFN_ELEMENTS(file_mod07)GT0THENBEGINhdf07_id=HDF_OPEN(file_mod07[0],/READ)water_vapor=MEAN(HDFSD_GETDATA(hdf07_id,'Water_Vapor',0,0),/NAN)HDF_CLOSE,hdf07_idENDIFELSEBEGINwater_vapor=5.0;热带默认值ENDELSEwater_vapor=MIN([MAX([water_vapor,4.0]),8.0]);观测天顶角（热带无海拔修正）file_mod03=FILE_SEARCH(input_path,'MYD03*'+STRMID(file_name,9,13)+'*.hdf')zenith_mean=(N_ELEMENTS(file_mod03)GT0)?MEAN(HDFSD_GETDATA(HDF_OPEN(file_mod03[0],/READ),'SensorZenith',0,0),/NAN):25.0;-------------------------步骤4：MODTRAN热带大气校正-------------------------modtran_pro,/INIT,$ATMOSPHERE='TROPICAL',$;热带大气模式ZENITH=zenith_mean,$WATER_VAPOR=water_vapor,$WAVELENGTH=[11.03],$/THERMAL,$OUTPUT_FILE=output_path+'modtran_'+file_name+'.dat'modtran_pro,/RUN,ERROR=mod_errIFmod_errNE0THENBEGINFPRINT,log_file,'错误：'+file_name+'-MODTRAN运行失败'CONTINUEENDIF;读取MODTRAN结果read_modtran,output_path+'modtran_'+file_name+'.dat',L_up,L_down,r_atm;热带高植被反射率修正rho_31=0.015;热带雨林热红外反射率Ls_31=L31-L_up-(L_down*rho_31)/(1-rho_31*r_atm+1e-8);-------------------------步骤5：热带比辐射率计算-------------------------;NDVI（高植被覆盖，NDVI>0.8）ndvi=(rho_b2-rho_b1)/(rho_b2+rho_b1+1e-8)ndvi[ndviLT-1ORndviGT1]=!VALUES.F_NAN;植被覆盖度（热带NDVI阈值）ndvi_min=0.1ndvi_max=0.8fv=(ndvi-ndvi_min)/(ndvi_max-ndvi_min+1e-8)fv[fvLT0]=0.0fv[fvGT1]=1.0fv[WHERE(FINITE(fv)EQ0)]=1.0;NaN设为全植被;热带雨林比辐射率eps_v=0.988;热带雨林植被eps_s=0.978;林下裸土d_eps=0.0015;混合修正项eps_31=fv*eps_v+(1-fv)*eps_s+d_eps;热带云掩膜（NDVI>0.8且亮温<290K判定为云）T31=14388/(11.03*ALOG(3.7418e16/(11.03^5*L31)+1))cloud_mask=(ndviGT0.8)AND(T31LT290)eps_31[cloud_mask]=!VALUES.F_NANLs_31[cloud_mask]=!VALUES.F_NAN;-------------------------步骤6：普朗克迭代（热带高温适配）-------------------------c1=3.7418e-16c2=1.4388e-2lambda31=11.03e-6Ts=FLTARR(SIZE(Ls_31,/DIMENSIONS))FORi=0,N_ELEMENTS(Ls_31)-1DOBEGINIF~FINITE(Ls_31[i])OR~FINITE(eps_31[i])ORLs_31[i]<=0THENBEGINTs[i]=!VALUES.F_NANCONTINUEENDIF;热带初始值（280~320K）T0=c2/(lambda31*ALOG(c1/(lambda31^5*L31[i])+1))T0=MIN([MAX([T0,280.0]),320.0])T1=T0delta=1.0iter_num=0WHILEABS(delta)GT0.01ANDiter_numLT50DOBEGINLb=c1/(lambda31^5*(EXP(c2/(lambda31*T1))-1))Ls_calc=eps_31[i]*LbdLdT=eps_31[i]*c1*c2/(lambda31^6*(EXP(c2/(lambda31*T1))-1)^2)*EXP(c2/(lambda31*T1))IFABS(dLdT)<1e-10THENBEGINdelta=0.0T1=!VALUES.F_NANBREAKENDIFdelta=(Ls_31[i]-Ls_calc)/dLdTdelta=MIN([MAX([delta,-1.0]),1.0])T1=T1+deltaiter_num=iter_num+1ENDWHILETs[i]=(iter_num>=50)?!VALUES.F_NAN:T1ENDFOR;-------------------------步骤7：结果保存-------------------------Ts_C=Ts-273.15Ts_C[Ts_CLT15ORTs_CGT40]=!VALUES.F_NAN;热带温度范围envi_map=ENVI_GET_MAP_INFO(FILE=file_mod02)ENVI_WRITE_ENVI_FILE,Ts_C,$FILENAME=output_path+file_name+'_LST.dat',$INTERLEAVE=0,$MAP_INFO=envi_map,$DESCRIPTION='热带森林LST（℃）_MODTRAN校正';-------------------------精度验证-------------------------IFFILE_TEST(input_path+'Tropical_meas_data.txt')THENBEGINmeas_data=READ_ASCII(input_path+'Tropical_meas_data.txt',DATA_COLUMNS=[3])meas_T=meas_data.DATA[0,*]lat=meas_data.DATA[1,*]lon=meas_data.DATA[2,*]ENVI_CONVERT_COORD,envi_map,lon,lat,x,y,/GEO_TO_IMAGEinv_T=FLTARR(N_ELEMENTS(x))FORj=0,N_ELEMENTS(x)-1DOBEGINix=FIX(x[j])iy=FIX(y[j])IFixGE0ANDixLTSIZE(Ts_C,1)ANDiyGE0ANDiyLTSIZE(Ts_C,2)THENinv_T[j]=Ts_C[iy,ix]ELSEinv_T[j]=!VALUES.F_NANENDFORvalid_idx=WHERE(FINITE(inv_T)ANDFINITE(meas_T),n_valid)IFn_validGT0THENBEGINrmse=SQRT(MEAN((inv_T[valid_idx]-meas_T[valid_idx])^2))FPRINT,log_file,file_name+'-RMSE：'+FORMAT='(F6.2)',rmse+'℃'ENDIFENDIFHEAP_FREE,data31,data32,L31,Ts,Ts_CPRINT,file_name+'处理完成'ENDFORFCLOSE,log_filePRINT,'热带森林LST反演完成！结果路径：',output_pathEND;复用前文read_modtran辅助函数PROread_modtran,modtran_file,L_up,L_down,r_atmL_up=0.0L_down=0.0r_atm=0.0IF~FILE_TEST(modtran_file)THENRETURNlun=GET_LUN()OPENR,lun,modtran_file,/GET_LINEline=''WHILE~EOF(lun)DOBEGINREADF,lun,lineIFSTRPOS(line,'UPWARDRADIANCE')NE-1THENBEGINREADF,lun,lineL_up=FLOAT(STRMID(line,20,15))*1e-6ENDIFIFSTRPOS(line,'DOWNWARDRADIANCE')NE-1THENBEGINREADF,lun,lineL_down=FLOAT(STRMID(line,20,15))*1e-6ENDIFIFSTRPOS(line,'ATMOSPHERICALBEDO')NE-1THENBEGINREADF,lun,liner_atm=FLOAT(STRMID(line,20,15))ENDIFENDWHILEFREE_LUN,lunEND二、中纬度农田MODISLST反演（完整IDL代码）场景适配：华北/美国中部农田、四季分明、混合像元为主idlPROMidLat_Farmland_LST;1.路径配置input_path='D:/MODIS_MidLat/Farmland/'output_path='D:/MODIS_MidLat/Result/'modtran_path='C:/MODTRAN6/'SETENV,'MODTRAN_DIR='+modtran_pathFILE_MKDIR,output_path;2.读取文件+日志files_mod02=FILE_SEARCH(input_path,'MOD021KM*.hdf')log_file=FOPEN(output_path+'Farmland_LST_Log.log',/WRITE)FPRINT,log_file,'中纬度农田LST反演日志|日期：'+SYSTIME();3.分季节参数配置season_params={summer:{atm:'MIDLAT_SUMMER',water:2.8,ndvi_min:0.1,ndvi_max:0.8,eps_v:0.987,eps_s:0.976,a3:1.35},winter:{atm:'MIDLAT_WINTER',water:1.8,ndvi_min:0.0,ndvi_max:0.5,eps_v:0.980,eps_s:0.973,a3:1.30}};=====================================批量处理=====================================FORf_idx=0,N_ELEMENTS(files_mod02)-1DOBEGINfile_mod02=files_mod02[f_idx]file_name=FILE_BASENAME(file_mod02,'.hdf')PRINT,'处理：',file_name;识别季节（从文件名提取月份：A2025001→第1天→1月，A2025180→第180天→6月）day_of_year=STRMID(file_name,9,3)month=CEIL(FLOAT(day_of_year)/30.5)IFmonthGE5ANDmonthLE9THENseason='summer'ELSEseason='winter'params=season_params.(season);-------------------------步骤1：数据读取+定标-------------------------hdf_id=HDF_OPEN(file_mod02,/READ,ERROR=err)IFerrNE0THENBEGINFPRINT,log_file,'错误：'+file_name+'-HDF打开失败'CONTINUEENDIFdata31=HDFSD_GETDATA(hdf_id,'EV_1KM_Emissive',0,10)data32=HDFSD_GETDATA(hdf_id,'EV_1KM_Emissive',0,11)gain31=HDFATTR_GET(hdf_id,'EV_1KM_Emissive','gain',10)offset31=HDFATTR_GET(hdf_id,'EV_1KM_Emissive','offset',10)data_b1=HDFSD_GETDATA(hdf_id,'EV_1KM_Reflective',0,0)data_b2=HDFSD_GETDATA(hdf_id,'EV_1KM_Reflective',0,1)HDF_CLOSE,hdf_idL31=data31*gain31+offset31rho_b1=data_b1*HDFATTR_GET(hdf_id,'EV_1KM_Reflective','gain',0)+HDFATTR_GET(hdf_id,'EV_1KM_Reflective','offset',0)rho_b2=data_b2*HDFATTR_GET(hdf_id,'EV_1KM_Reflective','gain',1)+HDFATTR_GET(hdf_id,'EV_1KM_Reflective','offset',1);-------------------------步骤2：MODTRAN大气校正（分季节）-------------------------modtran_pro,/INIT,$ATMOSPHERE=params.atm,$WATER_VAPOR=params.water,$ZENITH=28.0,$;中纬度平均天顶角WAVELENGTH=[11.03],$/THERMAL,$OUTPUT_FILE=output_path+'modtran_'+file_name+'.dat'modtran_pro,/RUN,ERROR=mod_errIFmod_errNE0THENCONTINUEread_modtran,output_path+'modtran_'+file_name+'.dat',L_up,L_down,r_atmrho_31=0.02;农田反射率Ls_31=L31-L_up-(L_down*rho_31)/(1-rho_31*r_atm+1e-8);-------------------------步骤3：分季节比辐射率-------------------------ndvi=(rho_b2-rho_b1)/(rho_b2+rho_b1+1e-8)ndvi[ndviLT-1ORndviGT1]=!VALUES.F_NANfv=(ndvi-params.ndvi_min)/(params.ndvi_max-params.ndvi_min+1e-8)fv[fvLT0]=0.0fv[fvGT1]=1.0eps_31=fv*params.eps_v+(1-fv)*params.eps_s+0.001;-------------------------步骤4：迭代反演-------------------------c1=3.7418e-16c2=1.4388e-2lambda31=11.03e-6Ts=FLTARR(SIZE(Ls_31,/DIMENSIONS))FORi=0,N_ELEMENTS(Ls_31)-1DOBEGINIF~FINITE(Ls_31[i])OR~FINITE(eps_31[i])THENBEGINTs[i]=!VALUES.F_NANCONTINUEENDIFT0=c2/(lambda31*ALOG(c1/(lambda31^5*L31[i])+1))T0=(seasonEQ'summer')?MIN([MAX([T0,280.0]),320.0]):MIN([MAX([T0,250.0]),290.0])T1=T0delta=1.0iter_num=0WHILEABS(delta)GT0.01ANDiter_numLT50DOBEGINLb=c1/(lambda31^5*(EXP(c2/(lambda31*T1))-1))Ls_calc=eps_31[i]*LbdLdT=eps_31[i]*c1*c2/(lambda31^6*(EXP(c2/(lambda31*T1))-1)^2)*EXP(c2/(lambda31*T1))delta=(Ls_31[i]-Ls_calc)/(dLdT+1e-10)delta=MIN([MAX([delta,-1.0]),1.0])T1=T1+deltaiter_num=iter_num+1ENDWHILETs[i]=(iter_num>=50)?!VALUES.F_NAN:T1ENDFOR;-------------------------结果保存-------------------------Ts_C=Ts-273.15Ts_C[(seasonEQ'summer')AND(Ts_CLT15ORTs_CGT35)]=!VALUES.F_NANTs_C[(seasonEQ'winter')AND(Ts_CLT-10ORTs_CGT20)]=!VALUES.F_NANenvi_map=ENVI_GET_MAP_INFO(FILE=file_mod02)ENVI_WRITE_ENVI_FILE,Ts_C,$FILENAME=output_path+file_name+'_LST.dat',$INTERLEAVE=0,$MAP_INFO=envi_mapFPRINT,log_file,file_name+'-季节：'+season+'-温度范围：'+STRTRIM(MIN(Ts_C,/NAN),2)+'~'+STRTRIM(MAX(Ts_C,/NAN),2)+'℃'HEAP_FREE,data31,L31,Ts,Ts_CENDFORFCLOSE,log_filePRINT,'中纬度农田LST反演完成！'END三、高频报错解决方案1.IDL读取HDF失败（报错：HDF_OPENerror/数据为空）原因&解决方案：报错现象原因解决方案HDF_OPEN:Erroropeningfile文件路径含中文/空格①将数据路径改为纯英文（如D:/MODIS_Data/，避免D:/MODIS数据/）；②文件名不含空格/特殊字符。HDFSD_GETDATA:Dataisempty波段索引错误核对MODIS波段索引：Band31=10、Band32=11（EV_1KM_Emissive），Band1=0、Band2=1（EV_1KM_Reflective）。HDFATTR_GET:Attributenotfoun