IDL中利用shape文件裁剪栅格图像

---

前言

本机环境：IDL 8.5 ENVI5.3

下载器要做数据预处理了，先裁剪原图像，方便后面处理。
主要利用了官方提供的方法

使用教程

利用shape文件裁剪栅格图像。
将完整代码键入IDL中，重置-编译-运行即可。

完整程序

pro RasterSubset
  COMPILE_OPT idl2
  ENVI,/restore_base_save_files
  envi_batch_init,LOG_FILE='batch.log'
 
  ;  多个文件夹的名字
  file_dic=['1','2','3','4']
  indexes = n_elements(file_dic) - 1
  for index=0,indexes  do begin
    ;  打开要裁剪的图像，分别是
    image_dir='F:\'+file_dic[index]+'\Landsat8' ;  图像所在文件夹路径
    shape_dir='F:\'+file_dic[index]+'\shape'  ;  shape文件文件夹路径
    roi_dir='F:\'+file_dic[index]+'\ROItest\' ;  生成的文件保存路径
    file_mkdir,roi_dir  ;  新建文件夹
    ;  在对应的目录下查找文件，数量为numfiles.
    ;  根据相应的文件格式修改过滤条件
    image_files=file_search(image_dir,'*blue.img',count=numfiles)  
    shapefile=file_search(shape_dir,'*.shp')
    if strlen(shapefile) eq 0 then return

    for i=0,numfiles-1 do begin
      image_file=image_files[i]
      print,image_file
      if strlen(image_file) eq 0 then return
      ENVI_OPEN_FILE, image_file, r_fid=fid, /no_interactive_query, /no_realize
      IF fid EQ -1 THEN  RETURN
      ;  生成保存文件的文件名和路径，根据自己的文件更改。
      dicnames=image_file.Split('\\')
      dicname=dicnames[-2].Split('_')
      out_name = roi_dir+dicname[-5]+'_'+dicname[-4]+'_'$
        +file_baseName(image_file,'.img')+'.img'    
      ;  调用官方的方法
      RasterSubsetViaShapefile,fid,shpFile=shapefile,outFile=out_name
      print, format='(%"-------%d: %d / %d -------")',index+1,i+1,numfiles
    endfor

  endfor
  ;tmp = DIALOG_MESSAGE('裁切结束!',/info)
  envi_batch_exit
END

;+
; :Description:
;    利用shapefile对栅格图像进行裁剪.
;
; :Syntax
;    RasterSubsetViaShapefile, Fid, shpFile=string, [pos=array],
;       [inside={0|1}], [outFile=｛string|variable}], [r_fid=variable]
;
; :Params:
;    Fid     -- 输入文件FID
;               注：可通过ENVI_OPEN_FILE、ENVI_SELECT、ENVIRasterToFID等获取
;
; :Keywords:
;    pos     -- 保留波段索引数组（可选），默认保留所有波段。
;    shpFile -- 用于裁剪的shapefile完整路径
;    inside  -- 保留shp文件外或内（可选，0或1），默认保留内部。
;               注：设置0时，保留外部；设置1时，保留内部。
;    outFile -- 裁减结果文件路径（可选）
;               注：如果不设置或设置为变量，则裁剪结果保存在临时目录中，outFile将保存输出文件名
;               注：如果设置为文件路径，则裁剪结果保存在指定路径中
;    r_fid   -- 返回裁剪结果文件FID，如果范围-1，则表示裁剪失败。
;
; :Author: duhj@esrichina.com.cn
;
; :Date: 2015年9月2日 15:33:38
;-
PRO RasterSubsetViaShapefile, Fid, shpFile=shpFile, pos=pos, $
  inside=inside, outFile=outFile, r_fid=r_fid

  COMPILE_OPT idl2
  ENVI, /RESTORE_BASE_SAVE_FILES
  ENVI_BATCH_INIT

  CATCH, err
  IF (err NE 0) THEN BEGIN
    CATCH, /CANCEL
    PRINT, 'ERROR: ' + !ERROR_STATE.MSG
    MESSAGE, /RESET
    RETURN
  ENDIF

  ENVI_FILE_QUERY, fid, ns=ns, nl=nl, nb=nb, $
    dims=dims, fname=fname, bnames=bnames

  IF ~N_ELEMENTS(pos)      THEN pos     = LINDGEN(nb)
  IF ~N_ELEMENTS(inside)   THEN inside  = 1
  IF ~KEYWORD_SET(outFile) THEN outFile = envi_get_tmp()

  ;读取shp文件的信息
  oshp=OBJ_NEW('IDLffShape',shpFile)
  IF ~OBJ_VALID(oshp) THEN RETURN
  oshp->GETPROPERTY,n_entities=n_ent,$ ;记录个数
    Attribute_info=attr_info,$ ;属性信息，结构体， name为属性名
    ATTRIBUTE_NAMES = attr_names, $
    n_attributes=n_attr,$ ;属性个数
    Entity_type=ent_type  ;记录类型

  iProj = ENVI_PROJ_CREATE(/geographic)
  ;自动读取prj文件获取投影坐标系
  potPos = STRPOS(shpFile,'.',/reverse_search)  ;
  prjfile = STRMID(shpFile,0,potPos[0])+'.prj'

  IF FILE_TEST(prjfile) THEN BEGIN
    OPENR, lun, prjFile, /GET_LUN
    strprj = ''
    READF, lun, strprj
    FREE_LUN, lun

    CASE STRMID(strprj, 0,6) OF
      'GEOGCS': BEGIN
        iProj = ENVI_PROJ_CREATE(PE_COORD_SYS_STR=strprj, $
          type = 1)
      END
      'PROJCS': BEGIN
        iProj = ENVI_PROJ_CREATE(PE_COORD_SYS_STR=strprj, $
          type = 42)
      END
    ENDCASE
  ENDIF

  oProj = ENVI_GET_PROJECTION(fid = fid)

  ;然后使用ROI进行掩膜统计
  roi_ids = !NULL
  FOR i = 0, n_ent-1 DO BEGIN
    ;
    ent = oshp->GETENTITY(i, /ATTRIBUTES) ;第i条记录

    ;如果ent不是多边形，则继续
    IF ent.SHAPE_TYPE NE 5 THEN CONTINUE

    N_VERTICES=ent.N_VERTICES ;顶点个数

    parts=*(ent.PARTS)

    verts=*(ent.VERTICES)
    ; 将顶点坐标转换为输入文件的地理坐标
    ENVI_CONVERT_PROJECTION_COORDINATES,  $
      verts[0,*], verts[1,*], iProj,    $
      oXmap, oYmap, oProj
    ; 转换为文件坐标
    ENVI_CONVERT_FILE_COORDINATES,fid,    $
      xFile,yFile,oXmap,oYmap

    IF (MIN(xFile) GE ns OR $
      MIN(yFile) GE nl OR $
      MAX(xFile) LE 0 OR $
      MAX(yFile) LE 0) AND i NE 0 THEN CONTINUE

    ;记录XY的区间，裁剪用
    IF i EQ 0 THEN BEGIN
      xmin = ROUND(MIN(xFile,max = xMax))
      yMin = ROUND(MIN(yFile,max = yMax))
    ENDIF ELSE BEGIN
      xmin = xMin < ROUND(MIN(xFile))
      xMax = xMax > ROUND(MAX(xFile))
      yMin = yMin < ROUND(MIN(yFile))
      yMax = yMax > ROUND(MAX(yFile))
    ENDELSE

    ;创建ROI
    N_Parts = N_ELEMENTS(Parts)

    FOR j=0, N_Parts-1 DO BEGIN
      roi_id = ENVI_CREATE_ROI(color=i,     $
        ns = ns ,  nl = nl)
      IF j EQ N_Parts-1 THEN BEGIN
        tmpFileX = xFile[Parts[j]:*]
        tmpFileY = yFile[Parts[j]:*]
      ENDIF ELSE BEGIN
        tmpFileX = xFile[Parts[j]:Parts[j+1]-1]
        tmpFileY = yFile[Parts[j]:Parts[j+1]-1]
      ENDELSE

      ENVI_DEFINE_ROI, roi_id, /polygon,    $
        xpts=REFORM(tmpFileX), ypts=REFORM(tmpFileY)

      ;如果有的ROI像元数为0，则不保存
      ENVI_GET_ROI_INFORMATION, roi_id, NPTS=npts
      IF npts EQ 0 THEN CONTINUE

      roi_ids = [roi_ids, roi_id]
    ENDFOR
  ENDFOR

  ;创建掩膜，裁剪后掩
  ENVI_MASK_DOIT,         $
    AND_OR = 2,           $
    IN_MEMORY=0,          $
    ROI_IDS= roi_ids,     $ ;ROI的ID
    ns = ns, nl = nl,     $
    inside=inside,        $ ;区域内或外
    r_fid = m_fid,        $
    out_name = envi_get_tmp()

  CASE inside OF
    0: out_dims=dims
    1: BEGIN
      ;进行掩膜 ---即不规则裁剪
      ;out_dims必须round，否则报错
      xMin = xMin >0
      xMax = ROUND(xMax) < (ns-1)
      yMin = yMin >0
      yMax = ROUND(yMax) < (nl-1)
      out_dims = [-1,xMin,xMax,yMin,yMax]
    END
  ENDCASE

  ENVI_MASK_APPLY_DOIT, FID = fid,      $
    POS = POS,                          $
    DIMS = out_dims,                    $
    M_FID = m_fid, M_POS = [0],         $
    VALUE = 0, r_fid = r_fid,           $
    OUT_NAME = outFile

  ; 掩膜文件ID移除
  ENVI_FILE_MNG, ID = m_fid,/REMOVE,/DELETE

  OBJ_DESTROY,oshp
END