PIOMAS月海冰厚度数据批量处理+建立地理信息（IDL的IGM建立地理信息）

下载的为PIOMAS海冰厚度的月数据，我感觉做了批量之后代码变复杂了一些，看不懂可以看之前的博客，读取单个文件代码比较简单。

首先我的程序分为两个部分，先把PIOMAS的数据提取出来，存储的格式为envi标准格式，将一年的数据存储在一起，即为一个数据12个波段。

pro read_piomas_monthly
nv = envi(/HEADLESS)
  infilepath = 'H:\mission\PIOMAS\monthly\'
  outfilename='H:\mission\PIOMAS\monthly\hdr\'
  CD,infilepath
  thesefiles = FILE_SEARCH('heff*',count = numFiles)
  ;解压
  FOR fidx=0,numFiles-1   DO BEGIN;
    s=strsplit(thesefiles[fidx],'.',/extract)
    if N_elements(s) eq 3 and  file_test (s[0]+'.'+s[1]) eq 0 then begin ;eq 0代表文件不存在
      FILE_GUNZIP, thesefiles[fidx], /VERBOSE
    endif
  ENDFOR
  nx=360 & ny=120;120 ;Arctic Parallel Ocean and sea Ice Model (POIM) grid
  mask_value=99999.0
  barwidth=200.0
  lon=fltarr(nx,ny) & lat=fltarr(nx,ny)
  data4=fltarr(nx,ny,12, /NOZERO);59
  kmt=intarr(nx,ny)
  thesefiles = FILE_SEARCH('heff*',count = numFiles)
  ;处理数据
  ;读取经纬度
  openr, lunin, 'H:\mission\LNG\PIOMAS\grid.dat',/Get_lun  ;打开数据文件
  readf, lunin, lon, format='(10f8.2)'
  readf, lunin, lat, format='(10f8.2)'
  close,lunin
  free_lun,lunin
  ;read grid mask (ocean levels)
  openr,lunin,'H:\mission\LNG\PIOMAS\io.dat_360_120.output',/Get_lun  ;打开数据文件
  readf, lunin, kmt,format='(360i2)'
  close,lunin
  free_lun,lunin
 
  
  FOR fidx=0,numFiles-1   DO BEGIN;
    s=strsplit(thesefiles[fidx],'.',/extract)
    if N_elements(s) eq 2 then begin ;eq 0代表文件不存在
      print,s
      ;读取数据
      if strcmp(s[1],'H2021') eq 1 then begin
        data4=fltarr(nx,ny,2, /NOZERO);只有两个月的数据
      endif else begin
        data4=fltarr(nx,ny,12, /NOZERO)
      endelse
      openr,lunin, thesefiles[fidx], /swap_if_big_endian,/Get_lun  ;打开数据文件'H:\mission\LNG\PIOMAS\monthly\heff.H2021'
      readu,lunin, data4
      close,lunin
      free_lun,lunin
      location=where(kmt lt 1,count)
      nband=size(data4)
      nband=nband[3]
      if count gt 0 then begin
        FOR band=0, nband-1 DO BEGIN
          temp=data4[*,*,band]
          temp[location]=0.0/0
          data4[*,*,band]=temp[*,*]
        ENDFOR
      endif
      print,'count',count
      ;存储无地理信息的数据
      ;存储头文件
      
      envi_setup_head,FNAME=outfilename+'piomas_monthly_'+s[1],NS=nx,NL=ny,NB=nband,DATA_TYPE=4,INTERLEAVE=0,/OPEN,/WRITE
      
      ;存储数据
      openw,lunou,outfilename+'piomas_monthly_'+s[1],/get_lun
      writeu,lunou,data4;
      close,lunou
      free_lun,lunou
      print,outfilename+'piomas_monthly_'+s[1],' finished'
    endif
  ENDFOR
nv.close
end

之后将标准格式转为带地理信息的tif，我使用的也是Georeference from IGM建立地理信息，但是这次换成代码的了。

首先说明我转了投影，我Georeference from IGM建立投影的详细信息，如下，以下的代码就为实现这个功能。

pro PIOMAS_geo_from_GLT
  infilepath = 'H:\mission\PIOMAS\monthly\hdr\'
  outfilename='H:\mission\PIOMAS\monthly\tif\'
  CD,infilepath
  thesefiles = FILE_SEARCH('piomas*',count = numFiles)
  ;打开、查询和读取Longitude数据
  envi_open_file, 'H:\mission\PIOMAS\lon', r_fid = x_fid
  envi_file_query, x_fid, dims = x_dims
  print, x_dims
  x_pos = 0
  data = envi_get_data(fid = x_fid, dims = x_dims, pos = x_pos)

  ;打开、查询和读取Latitude数据
  envi_open_file, 'H:\mission\PIOMAS\lat', r_fid = y_fid
  envi_file_query, y_fid, dims = y_dims
  print, y_dims
  y_pos = 0
  data = envi_get_data(fid = y_fid, dims = x_dims, pos = y_pos)
  
  ;建立GLT
  ;Make the GLT
  outname = outfilename  +'PIOMAS_glt'
  pixel_size = 25000d;12500d;250.0d;
  rotation = 0.0

  i_proj = envi_proj_create(/geographic, datum="WGS-84")

  params = [ 0,0,90,35]  ;PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Latitude_Of_Origin",90.0],
  datum = 'WGS-84'
  name = 'PIOMAS azimuthal equal area'
  units = envi_translate_projection_units('Meters')

  o_proj = envi_proj_create(type=11, name = name, datum = datum, params = params, units = units) ;11;36 ;31-Stereographic_North_Pole
  print, o_proj

  envi_glt_doit, i_proj=i_proj, $

    o_proj=o_proj, out_name = outname, dims = dims, $

    pixel_size = pixel_size, r_fid = glt_fid, $

    rotation = rotation, x_fid = x_fid, y_fid = y_fid, $

    x_pos = x_pos, y_pos = y_pos
  if glt_fid eq -1 then begin
    print,'GLT建立失败'
  endif
  ;利用GLT建立地理信息

  if glt_fid eq -1 then return
  
  FOR fidx=0,numFiles-1   DO BEGIN;
    s=strsplit(thesefiles[fidx],'.',/extract)
    if n_elements(s) gt 1 then begin
      continue
    endif
    ;打开、查询和读取image数据
    envi_open_file, thesefiles[fidx], r_fid = fid
    envi_file_query, fid, dims = dims,NB=NB
    PRINT, dims
    ;  dims[1] = 24
      print, thesefiles[fidx],NB
    pos = lindgen(nb)
    data = envi_get_data(fid = fid, dims = dims, pos = pos)
    
    envi_file_query, glt_fid, dims = glt_dims

    out_name =  outfilename+thesefiles[fidx]+ '_georef'

    envi_doit, 'envi_georef_from_glt_doit', fid = fid, $

    glt_fid = glt_fid, glt_dims = glt_dims, out_name = out_name, pos = pos, $

    r_fid = r_fid,BACKGROUND=-1
    
    Tif_File = outfilename+thesefiles[fidx]+ '_georef_tif.tif' ;输出的文件名
    ENVI_OUTPUT_TO_EXTERNAL_FORMAT,FID=r_fid,Dims=glt_dims,POS=pos,Out_Name=Tif_File,/TIFF

  ENDFOR
  

  

end