IDL矢量/掩膜裁剪影像(ENVI API编程)

一、 矢量裁剪

1、问题描述

当我们想要用一个矢量文件去裁剪影像时，可以使用ENVIVectorMaskRaster函数。

maskedRaster = ENVIVectorMaskRaster(raster, shp_file)

图中标记红圈处是我们裁剪出的影像，非常小的一块。但是却发现用ENVI-IDL裁剪出的影像仍然保留着原来影像的行列数，并不像ENVI中Subset Data from ROIs工具一样可以自行适应较小矢量文件的行列数。
那么这个时候就需要我们重新定义影像的空间范围。

2、重定义空间范围

2.1定义范围

Grid = ENVIGridDefinition(CoordSys, $
  PIXEL_SIZE=[9.186D,9.186D], $
  TIE_POINT_PIXEL=[0.0D,0.0D], $
  TIE_POINT_MAP=[3075299.7946D,1246937.9905D], $
  NROWS=Raster.NROWS, $
  NCOLUMNS=Raster.NCOLUMNS)

CoordSys:栅格或矢量的坐标系信息。

PIXEL_SIZE：像元大小。

TIE_POINT_PIXEL：Specify a two-element array with the map coordinates of the TIE_POINT_PIXEL location, as follows:[xmin, ymax]。
即指定MAP坐标的xmin和ymax
If you set this property, you must also specify NROWS, NCOLUMNS, and PIXEL_SIZE.
如果设定了[xmin, ymax]则还要指定影像的行列数，像元大小。

NROWS：影像行数。

NCOLUMNS：影像列数。

2.2转换到到新的空间范围

reGrid = ENVISpatialGridRaster(maskedRaster,GRID_DEFINITION=Grid)

3、矢量裁剪实例

上面提到了定义新的空间范围需要裁剪后的影像行列号，结合实例来说明行列号的获取。

pro ENVIVectorMaskRaster_
  compile_opt idl2, hidden
  
  file = 'E:\FY3D_process\FY3D_MERSI_data\MERSI2_20190403_0445_B3_4TOA.dat'
  file2 = 'F:\ENVI_Tempfn\Shp'
  
  raster = e.OpenRaster(file)
  file_shp = e.OpenVector(file2)
  
  maskedRaster = ENVIVectorMaskRaster(raster, file_shp)
  outFile = e.GetTemporaryFilename()
  maskedRaster.Export, outFile, 'ENVI', DATA_IGNORE_VALUE=0
  ;这里必须要构建一个临时文件，设置忽略值为0然后输出一下，否则裁剪后是矢量真实边界的外接矩形。
  
  x = (file_shp.data_range[3]-file_shp.data_range[1])/250;行数
  y = (file_shp.data_range[2]-file_shp.data_range[0])/250 ;列数
  
  Grid = ENVIGridDefinition(file_shp.Coord_Sys, $
  PIXEL_SIZE=[250.0D,250.0D], $  ;像元大小250M
  TIE_POINT_PIXEL=[0.0D,0.0D], $
  TIE_POINT_MAP=[file_shp.data_range[0],file_shp.data_range[3]], $
  NROWS=x, $
  NCOLUMNS=y)
  
  outFile = 'F:\ENVI_Tempfn\sub'
  regrid = ENVISpatialGridRaster(maskedRaster,GRID_DEFINITION=Grid)
  regrid.export,outFile,'envi'
  
  end

file_shp是打开的矢量文件，file_shp.data_range可以获取影像的左下，右上的Map坐标(x,y)，单位是M。那么影像的长和宽就可以由以下计算得出(file_shp.data_range[3]-file_shp.data_range[1]),(file_shp.data_range[2]-file_shp.data_range[0]).然后再除以像元的大小，就得到了影像的行列数。

试验一下不规则shp：

满足要求。

二、 掩膜裁剪

1、问题描述

我们进行去云处理时，通常会根据云在某些波段TOA的特征来设定阈值，然后剔除。剔除操作可以通过设定一个二进制掩膜去裁剪原影像。
那么步骤就分为两步：创建掩膜和裁剪影像。

2、创建掩膜

raster = e.OpenRaster(file)
subset = ENVISubsetRaster(raster, BAND=[0])
;从读取的影像中摘出你要作为依据设定阈值的波段
threshold = [0.3D]
binaryMask = ENVIBinaryLTThresholdRaster(raster threshold); 大于0.3的将被遮住

threshold = [0.3D]
binaryMask = ENVIBinaryGTThresholdRaster(raster, threshold); 小于0.3的将被遮住

ranges = [[0.3d,1]]
maskedRaster = ENVIDataValuesMaskRaster(raster, ranges) ;介于0.30-1之间的将被遮住

同样也可以这样创建掩膜，但上面三种创建的形式优势在于可以针对每一个波段设置阈值，如
threshold = [0.3D,0.5D] 官方帮助

mask = (raster.GetData(BAND=0) GE 0.3)

3、裁剪影像

After_maskRaster = ENVIMaskRaster(raster, binaryMask)

4、掩膜裁剪实例

这里我读取了MERSI2的band3和4的TOA。并设定band3的TOA大于0.3为云。

pro cloudmask
  compile_opt idl2
  e = ENVI()
  tic

  file = 'E:\FY3D_process\FY3D_MERSI_data\MERSI2_20190403_0445_B3_4TOA.dat'
  raster = e.OpenRaster(file)
  
  subset = ENVISubsetRaster(raster,Band=[0])
  threshold = [30.]
  binaryMask = ENVIBinaryLTThresholdRaster(subset, threshold)
  
  ;outFile = e.GetTemporaryFilename()
  ;binaryMask.Export, outFile, 'ENVI'
  
  outfile = 'F:\ENVI_Tempfn\cloud.dat'
  After_maskRaster = ENVIMASKRaster(raster, binaryMask)
  After_maskRaster.Export, outfile, 'ENVI', DATA_IGNORE_VALUE = 0 

  toc
end

去云效果图如下（加了一个灰度分割）：

统计一下掩膜后影像的值：

  Band3 = raster.GetData(BANDS=0, PIXEL_STATE=pixelState)
  
  pos = WHERE(pixelState EQ 0, count)
  IF (count NE N_Elements(Band3)) THEN Band3 = Band3[pos]
  
  Print, 'Band3最小值:', Min(Band3)
  Print, 'Band3最大值:', Max(Band3)
  
  Band3最小值:       3.0893782
  Band3最大值:       29.999999

参考资料：官方帮助MASK