【IDL代码库】如何用IDL获取栅格范围平均海拔

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