利用IDL批量将栅格数据从地理坐标系投影到投影坐标系

最新推荐文章于 2023-10-22 20:45:19 发布

Rosers™

最新推荐文章于 2023-10-22 20:45:19 发布

阅读量854

点赞数

文章标签：经验分享

本文链接：https://blog.csdn.net/weixin_44308329/article/details/125278708

版权

很多遥感影像下载后是地理坐标系的，但是在很多计算中，是需要投影才可以进行，在这里以NASADEM投影为例

将一副已经投影好的影像作为参考，引入投影参数

  ;定义投影
  File = 'C:\Users\DELL\Desktop\test\ty\DEM_Mercator.tif'
  ENVI_OPEN_FILE, File ,r_fid=fid
  o_proj = ENVI_GET_PROJECTION(FID = fid)
  ; o_proj = ENVI_PROJ_CREATE(/geographic)

批量读入数据

 CD,'C:\Users\DELL\Desktop\test\'
  BINFiles = FILE_SEARCH("*.tif")
  FileCount = N_ELEMENTS(BINFiles)
  IF FileCount EQ 0 THEN RETURN

逐景影像投影

 FOR NX = 0,FileCount -1 DO BEGIN
    FileName = BINFiles[NX]
    ENVI_OPEN_FILE, FileName ,r_fid=fid
    IF (fid EQ -1)THEN BEGIN
      ENVI_BATCH_EXIT
      RETURN
    ENDIF
    ; Setupthe values for the keywords
    ENVI_FILE_QUERY, fid[0], dims=dims, nb=nb
    pos = LINDGEN(nb)
;    indexstr = STRPOS(FileName,'.tif')
;    out_name = STRMID(FileName,0,indexstr) +'_Geo.tif'
    out_name = 'C:\Users\DELL\Desktop\test\result\'+STRMID(FileName,0,7) +'_Geo.tif'
    ;设置输出像元大小
    o_pixel_size = [30,30];


    ENVI_CONVERT_FILE_MAP_PROJECTION, fid=fid, $
      pos=pos, dims=dims, o_proj=o_proj, $
      o_pixel_size=o_pixel_size, $
      out_name=out_name, warp_method=2, $
      resampling=0, background=-9999
    ENVI_FILE_MNG,id = fid,/remove
    
    print,FileName

完整代码：

PRO CONVERT_FILE_MAP_PROJECTION
  ; Firstrestore all the base save files.
  COMPILE_OPT IDL2
  ENVI,/restore_base_save_files
  e = ENVI(/headless)
  ;Initialize ENVI and send all errors
  ; andwarnings to the file batch.txt
  ENVI_BATCH_INIT, log_file='batch.txt'
  ; Open theinput file
  ;定义投影
  File = 'C:\Users\DELL\Desktop\test\ty\DEM_Mercator.tif'
  ENVI_OPEN_FILE, File ,r_fid=fid
  o_proj = ENVI_GET_PROJECTION(FID = fid)
  ; o_proj = ENVI_PROJ_CREATE(/geographic)

  CD,'C:\Users\DELL\Desktop\test\'
  BINFiles = FILE_SEARCH("*.tif")
  FileCount = N_ELEMENTS(BINFiles)
  IF FileCount EQ 0 THEN RETURN

  FOR NX = 0,FileCount -1 DO BEGIN
    FileName = BINFiles[NX]
    ENVI_OPEN_FILE, FileName ,r_fid=fid
    IF (fid EQ -1)THEN BEGIN
      ENVI_BATCH_EXIT
      RETURN
    ENDIF
    ; Setupthe values for the keywords
    ENVI_FILE_QUERY, fid[0], dims=dims, nb=nb
    pos = LINDGEN(nb)
;    indexstr = STRPOS(FileName,'.tif')
;    out_name = STRMID(FileName,0,indexstr) +'_Geo.tif'
    out_name = 'C:\Users\DELL\Desktop\test\result\'+STRMID(FileName,0,7) +'_Geo.tif'
    ;设置输出像元大小
    o_pixel_size = [30,30];


    ENVI_CONVERT_FILE_MAP_PROJECTION, fid=fid, $
      pos=pos, dims=dims, o_proj=o_proj, $
      o_pixel_size=o_pixel_size, $
      out_name=out_name, warp_method=2, $
      resampling=0, background=-9999
    ENVI_FILE_MNG,id = fid,/remove
    
    print,FileName
  ENDFOR
  ; ExitENVI
  ENVI_BATCH_EXIT
  print,'finish'
END