FLAASH已提供IDL接口(http://blog.sina.com.cn/s/blog_764b1e9d0102xxrk.html),但在FLAASH大气校正中需要输入研究区平均海拔,我们已经介绍过如何利用ENVI的统计功能获取研究区高程(http://blog.sina.com.cn/s/blog_764b1e9d0101drqg.html),下面介绍如何一种使用IDL获取高程方法,以下为功能源码及测试代码:
;+
; :description:
; :获取栅格范围的平均海拔,单位为Km
; 失败则返回-1,否则返回浮点型栅格范围平均海拔
; :author: Hanzt
; :2017-11-28
; -
function Get_Raster_Msl,raster
;编译规则优化
compile_opt idl2
;获取envi对象
if isa(e,'envi') then e=envi(/current) else e=envi()
;获取ENVI自带~900m分辨率DEM数据绝对路径(如果有更高分辨率DEM也可替换)
root_Dir=e.root_dir
demFile=filepath('GMTED2010.jp2',$
root_Dir=root_Dir,subdirectory=['data'])
;如果DEM不存在,则手动选DEM数据
if ~file_test(demfile) then begin
title='Select a DEM file which covered the extent of the input raster.'
demFile=envi_pickfile(title=title)
if demFile eq '' then return,-1
endif
;打开DEM,并记录错误信息err0
demRaster=e.openraster(demFile,error=err0)
;获取栅格行列数
spatialref1=raster.spatialref
ns=raster.ncolumns
nl=raster.nrows
;获取栅格四个角点投影坐标,mapx,mapy
spatialref1.convertfiletomap,[0,ns,ns,0],[0,0,nl,nl],mapx,mapy
coordSys1=envicoordsys(coord_sys_str=spatialref1.coord_sys_str)
;将mapx,mapy转换到DEM栅格投影坐标,用于DEM裁剪
spatialref2=demraster.spatialref
coordSys2=envicoordsys(coord_sys_str=spatialref2.coord_sys_str)
coordSys1.convertmaptomap,mapx,mapy, Map2X, Map2Y, coordSys2
;获取范围最值
x_max=max(Map2X,min=x_min)
y_max=max(Map2Y,min=y_min)
;获取裁剪范围
sub_rect=[x_min,y_min,x_max,y_max]
;用栅格范裁剪DEM数据
rasterSub=envisubsetraster( $
demRaster, $
spatialref=spatialref2, $
sub_rect=sub_rect, $
error=err1)
;如果裁剪失败,则返回-1
if rasterSub eq !null then begin
msg='No intersections between the input raster and the dem raster.'
!null=dialog_message(msg,/info)
return,-1
endif
;裁剪后DEM统计
Statistics =ENVIRasterStatistics(rasterSub,error=err2)
;如果失败返回-1,否则返回平高程
if err0 ne '' || err1 ne '' || err2 ne '' then $
return,-1 else $
return,float(Statistics['MEAN']/1000)
end
;测试代码
pro test_Get_Raster_Msl
compile_opt idl2
if isa(e,'envi') then e=envi(/current) else e=envi()
file='C:\Program Files\Exelis\ENVI53\data\qb_boulder_msi'
raster=e.openraster(file)
msl=get_raster_msl(raster)
;如果成功则控制台打印
if msl.typecode eq 4 then print,'Mean Sea Leavel:',msl
end