IDL 对CMIP6数据进行投影+地理信息

 这一步是基于IDL读取CMIP6的NC数据的下一步，之前那一步生成了经度和纬度及变量三个矩阵tif文件。

以下新版本更新了释放内存，在每次年循环后面增加释放内存，以防内存不足，以及处理速度越来越慢。

pro awi_geo_from_GLT
compile_opt idl2
  nv = envi(/HEADLESS)
for yr=2017,2020 do begin
  infilepath ='H:\mission\LNG\Model\MPI-ESM1-2-LR\ssp245\daily\sicon\hdr\'+STRING(yr, FORMAT='(I04)') 
  outfilename='H:\mission\LNG\Model\MPI-ESM1-2-LR\ssp245\daily\sicon\area\';
  CD,infilepath
  thesefiles = FILE_SEARCH('si*',count = numFiles);siconc sithick需要替换
  print,'一共文件： ',numFiles
  ;打开、查询和读取Longitude数据
  envi_open_file,'H:\mission\LNG\Model\MPI-ESM1-2-LR\ssp245\daily\sicon\lon' , r_fid = x_fid ;'H:\mission\LNG\Model\AWI-CM-1-1-MR\ssp370\scon\daily\lon'
  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\LNG\Model\MPI-ESM1-2-LR\ssp245\daily\sicon\lat', r_fid = y_fid; 'H:\mission\LNG\Model\AWI-CM-1-1-MR\ssp370\scon\daily\lat'
  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  +'awi_glt'
  pixel_size = 25000d;12500d;250.0d;
  rotation = 0.0

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

  params = [0,0,70,-45];后面是之前和AMSR2配准建立的，不过感觉不太准[ 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 = '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-azimuthal equal area;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
    str=STRSPLIT(thesefiles[fidx], '_', /EXTRACT)
    print, thesefiles[fidx],NB,STRMID(str[-1],0,4);STRMID(thesefiles[fidx],42,4)
    pos = lindgen(nb)
    data = envi_get_data(fid = fid, dims = dims, pos = pos)

    envi_file_query, glt_fid, dims = glt_dims

    out_name =  outfilename+STRMID(str[-1],0,4)+'\'+thesefiles[fidx]+ '_georef';
    print,out_name

    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+STRMID(str[-1],0,4)+'\'+thesefiles[fidx]+ '_georef_tif.tif' ;输出的文件名 
    ENVI_OUTPUT_TO_EXTERNAL_FORMAT,FID=r_fid,Dims=glt_dims,POS=pos,Out_Name=Tif_File,/TIFF

  ENDFOR
  ;释放内存

;释放内存
 
;获取当前内存中的所有文件的fid
 
  fids = envi_get_file_ids()
 
  ;获取数组的大小
 
  size = size(fids)
 
  length = size[1]
 
  ;循环释放内存中的文件
 
  ;主要就是envi_file_mng这个函数，其中id是打开文件的id，另外还有两个关键字
 
  ;remove是只从内存中移除，​delete不仅从内存中移除也从硬盘上删除，大家要              ;慎用
 
  for i = 0L, length-1 do begin
 
       envi_file_mng,id = fids[i],/remove
 
  endfor
 
endfor
nv.close
end