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

以你之姓@ 2022-11-13 00:53 249阅读 0赞

下载的为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建立投影的详细信息，如下，以下的代码就为实现这个功能。

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L1hpbmVtYUNoZW4_size_16_color_FFFFFF_t_70][]![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L1hpbmVtYUNoZW4_size_16_color_FFFFFF_t_70 1][]

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

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L1hpbmVtYUNoZW4_size_16_color_FFFFFF_t_70]: https://img-blog.csdnimg.cn/20210325095551389.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L1hpbmVtYUNoZW4=,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L1hpbmVtYUNoZW4_size_16_color_FFFFFF_t_70 1]: /images/20221022/7de883bacd814f7a85c6a2115b450222.png