后台很多朋友私信我,对如何进行均值替代反演无缺失值的LST产品感兴趣,手动在ENVI中处理太麻烦了,所以给出批量利用C2数据集生产LST的IDL代码,以供参考。
对原理实现不了解的朋友可以看我的上一篇博文,Landsat 计算LST(地表温度)——没有大气剖面参数计算器怎么办
代码在IDL8.8版本编译通过,以下给出代码实现:
pro cal_lst_from_c2,file_path=file_path,output_path=output_path,e=e
; 本pro实现根据Landsat C2数据集ST中间过程数据,采用均值替代法生产无缺失值的LST产品
; file_path : Landsat C2数据集路径,string型
; output_path: 输出文件的绝对路径,string型
; e : ENVI
COMPILE_OPT idl2
e=envi(/headless)
; 按自己的需求替换输入输出文件路径
file_path = 'E:\TEST\LC08_L2SP_127033_20180725_20200831_02_T1\'
output_path = 'E:\TEST\output\127033_20180725_LST.dat'
;计算上下行辐射、辐照强度
up_file = file_search(file_path,'*_URAD.tif')
down_file = file_search(file_path,'*_DRAD.tif')
tirs_file = file_search(file_path,'*_TRAD.tif')
upTask = enviTask('PixelwiseBandMathRaster')
upTask.input_raster = e.openraster(up_file[0])
upTask.DATA_IGNORE_VALUE = -9999
upTask.EXPRESSION = 'b1*0.001'
upTask.execute
downTask = enviTask('PixelwiseBandMathRaster')
downTask.input_raster = e.openraster(down_file[0])
downTask.DATA_IGNORE_VALUE = -9999
downTask.EXPRESSION = 'b1*0.001'
downTask.execute
tirsTask = enviTask('PixelwiseBandMathRaster')
tirsTask.input_raster = e.openraster(tirs_file[0])
tirsTask.DATA_IGNORE_VALUE = -9999
tirsTask.EXPRESSION = 'b1*0.001'
tirsTask.execute
;计算透过率、比辐射率、比辐射率偏差
atran_file = file_search(file_path,'*_ATRAN.tif')
emis_file = file_search(file_path,'*_EMIS.tif')
emsd_file = file_search(file_path,'*_EMSD.tif')
atranTask = enviTask('PixelwiseBandMathRaster')
atranTask.input_raster = e.openraster(atran_file[0])
atranTask.DATA_IGNORE_VALUE = -9999
atranTask.EXPRESSION = 'b1*0.0001'
atranTask.execute
emisTask = enviTask('PixelwiseBandMathRaster')
emisTask.input_raster = e.openraster(emis_file[0])
emisTask.DATA_IGNORE_VALUE = -9999
emisTask.EXPRESSION = 'b1*0.0001'
emisTask.execute
emsdTask = enviTask('PixelwiseBandMathRaster')
emsdTask.input_raster = e.openraster(emsd_file[0])
emsdTask.DATA_IGNORE_VALUE = -9999
emsdTask.EXPRESSION = 'b1*0.0001'
emsdTask.execute
;波段合成
stackTask = envitask('BuildLayerStack')
stackTask.input_rasters =[tirsTask.output_raster,upTask.output_raster,atranTask.output_raster,downTask.output_raster]
stackTask.execute
;统计比辐射率均值
emis_mean = (ENVIRasterStatistics(emisTask.output_raster)).mean
emsd_mean = (ENVIRasterStatistics(emsdTask.output_raster)).mean
total_emis = emis_mean[0]+emsd_mean[0]
;计算黑体辐射亮度
btTask = enviTask('PixelwiseBandMathRaster')
btTask.input_raster = stackTask.output_raster
btTask.DATA_IGNORE_VALUE = -9999
btTask.EXPRESSION = '(b1-b2-b3*(1-'+strcompress(string(total_emis))+')*b4)/(b3*'+strcompress(string(total_emis))+')'
btTask.execute
;反演地表温度
lstTask = enviTask('PixelwiseBandMathRaster')
lstTask.input_raster = btTask.output_raster
lstTask.DATA_IGNORE_VALUE = -9999
lstTask.EXPRESSION = '(1321.08)/alog(774.89/b1+1)-273.15'
lstTask.OUTPUT_RASTER_URI = output_path
lstTask.execute
end
欢迎留言讨论。