【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