IDL实现经纬度转行列号影像裁切

最近在做毕业设计，其中有一步需要利用经纬度裁切区域，作为IDL小白，看过好几个博客，实验了很多次，终于成功了，特此写下人生中第一篇csdn博客，纪念一下。特别感谢一下我的一位同学，他的博客写的很出色，IDL学的也特别好，毕设期间经常麻烦他。看到的小伙伴可以去看看他的博客，干货很多。https://www.ixxin.cn/

我用的IDL版本是IDL85，ENVI5.3 SP1

不过这个代码应该是各种版本基本都可以，对于FLAASH接口，就得用5.3以上才可以了

第一步是经纬度转行列号

PRO band_subset
  COMPILE_OPT IDL2
 
  e = envi(/headless)
 
  file = ‘F:\2013_04\OutDir\MYD02QKM.A20130403_radiance_georef.dat’
  dir = ‘F:\2013_04\out’
  name = file_basename(file)
  oFile = dir+’\’+name+’_subset.dat’   ;设置输出文件的名称
 
  MULRaster = envi.OpenRaster(File)
 
  Lon=[119,121]
  Lat=[36.5,35]
 
  MapXs = fltarr(2)
  MapYs = fltarr(2)
  FileXs = fltarr(2)
  FileYs = fltarr(2)
 
 ;经纬度转行列号
  for i = 0 ,1 do begin
    spatialRef1.ConvertLonLatToMap, Lon[i], Lat[i], MapX, MapY
    spatialRef1.ConvertMapToFile, MapX, MapY, FileX, FileY
    FileXs[i] = FileX
    FileYs[i] = FileY
    MapXs[i] = MapX
    MapYs[i] = MapY

  endfor

 ;由于我最终目的是利用行列号规则裁切，所以要使用round( )函数取整

 ;上面提到的那个同学给我讲解过round( )和fix( )取整机制的区别，小伙伴可以自行在IDL中试一试

  FileXs = round(FileXs)
  FileYs = round(FileYs)
 
  FileX1 = FileXs[0]
  FileX2 = FileXs[1]
 
  FileY1 = FileYs[0]

  FileY2 = FileYs[1]

到这一步，经纬度已经转为行列号了，下面是裁切部分，要用到ENVI函数库的一个函数

ENVISubsetRaster(Raster [, Keywords=value])（下面是ENVI帮助文档的函数介绍）

;这一步获取裁切区域

Subset = ENVISubsetRaster(MULRaster,SUB_RECT=[FileX1, FileY1, FileX2, FileY2])

；输出裁切后的影像数据

e.ExportRaster,Subset, oFile,’dat’

【注】ENVISubsetRaster（ ）的keywords里面包含bands，例如：

Subset = ENVISubsetRaster(MULRaster,SUB_RECT=[FileX1, FileY1, FileX2, FileY2]，bands = 0)下面的裁切过程将只对第一个波段进行，如果不指定这个参数，默认会对所有可用波段进行裁切，想要深入了解，建议去查一下这个函数的帮助文档，里面有各种example以及详细的说明。

这是一个完整的程序，希望能够帮助到有需要的小伙伴，以后有新的进步会继续更博客哒*^_^*

转载自：https://blog.csdn.net/Artemis_R/article/details/80159787