IDL CMIP6 NC格式数据处理

小灰灰 2022-10-08 02:28 297阅读 0赞

下载了一些模型包括以下：

<table style="width:228pt;"> 
 <tbody> 
  <tr> 
   <td style="vertical-align:middle;width:228pt;"><span style="color:#000000;">AWI-CM-1-1-MR</span></td> 
  </tr> 
  <tr> 
   <td style="vertical-align:middle;"><span style="color:#000000;">CMCC.CMCC-CM2-HR4</span></td> 
  </tr> 
  <tr> 
   <td style="vertical-align:middle;"><span style="color:#000000;">CNRM-CERFACS.CNRM-CM6-1-HR</span></td> 
  </tr> 
  <tr> 
   <td style="vertical-align:middle;"><span style="color:#000000;">NOAA-GFDL.GFDL-CM4</span></td> 
  </tr> 
  <tr> 
   <td style="vertical-align:middle;"><span style="color:#000000;">MOHC.HadGEM3-GC31-MM</span></td> 
  </tr> 
 </tbody> 
</table>

下载的是siconc和sithick两个变量，而且存在SIMon，和SIday两种frequency的。

所以要求程序能够自动处理不同模型、变量不同、频率不同的情况。

程序只把数据提取成envi标准格式，并且也存了envi的hdr，文件名需要保留除日期之外的全部，再加上日期/月份，存储的文件是月/日的。

最后也保存了lon和lat

我觉得比较难做的主要是，各种日期的识别，因为有考虑闰年的，不考虑闰年的之类，因此都单独用了函数来计算。

pro read_CMIP6_model_data
      infilepath = 'H:\mission\LNG\Model\GFDL-CM4\ssp585\daily\sicon'
      outfilepath = 'H:\mission\LNG\Model\GFDL-CM4\ssp585\daily\sicon\hdr\'
      nv = envi(/HEADLESS)
      ;  geo=read_tiff('H:\mission\LNG\mask\smos_geo2.tif',geotiff=geoinfo);
      ;  mask=read_tiff('H:\mission\LNG\mask\smos_landmask.tif');
      CD,infilepath
      thesefiles = FILE_SEARCH('*2075*.nc',count = numFiles)
      print,'number of files found:',numFiles
      ;处理数据
    ;  avresult=intarr(165,12+1) ;1850-2014;12行：12个月
    ;  avresult[*,*]=-1 ;缺数据都为0
    ;  avresult[*,0]=indgen(165)+1850 ;附上年份
      FOR fidx=0, numFiles-1 DO BEGIN; ;不要忘记修改这个
        filename=file_basename(thesefiles[fidx],'.nc')
        str=STRSPLIT(filename, '_', /EXTRACT) ;siconc_SIday_AWI-CM-1-1-MR_historical_r1i1p1f1_gn_18500101-18501231
        ;获取的是一期的数据
        file_ID = NCDF_OPEN(theseFiles[fidx],/NOWRITE)   ;open netCDF file for READ only
        Tag = NCDF_INQUIRE(file_ID)
        timeid = NCDF_VARID(file_ID,'time')             ;initial_time0initial_time0_hours  initial_time0
        NCDF_VARGET, file_ID, timeid, time              ;获取 时间维变量
        nt = N_ELEMENTS(time)
        NCDF_ATTGET,file_ID, timeid,'calendar',timeinfo
        timeinfo=string(timeinfo)
        NCDF_ATTGET,file_ID, timeid,'units',timeunits
        timeunits=string(timeunits)
        timeunits=STRSPLIT(timeunits, ' ', /EXTRACT)
        
        ;print,time
        latid = NCDF_VARID(file_ID,'lat')               ;g0_lat_2
        if latid eq -1 then begin
          latid = NCDF_VARID(file_ID,'latitude')  
        endif
        NCDF_VARGET, file_ID, latid, latitude           ;获取 时间维变量
        nlat = N_ELEMENTS(latitude)                     ;纬向维长度
        ;print,'dsklfjdsjfkldjsfkljl'+latitude
        lonid = NCDF_VARID(file_ID,'lon') 
        if lonid eq -1 then begin
          lonid = NCDF_VARID(file_ID,'longitude')
        endif               ;g0_lon_3
        NCDF_VARGET, file_ID, lonid, longitude           ;获取 时间维变量
        nlon = N_ELEMENTS(longitude)                     ;径向维长度
    
    
        var_ID = NCDF_VARID(file_ID,str[0] )  ;获取变量IDvar_result.Name  siconc'sithick'
        NCDF_VARGET, file_ID, Var_ID, var_Data        ;获取变量数据
        
        NCDF_ATTGET,file_ID, var_ID,'missing_value',varinfo
        varinfo=float(varinfo)
        
        ;大于等于0为海冰有效值，缺失值为-0.2，陆地为-0.1
    
        ;    var_Data=reverse(var_Data,2)
        print,max(var_Data)
    ;    ys=STRMID(str[6],0,8);20150101-20201231
    ;    ye=STRMID(str[6],9,8)
        timeunits=STRSPLIT(timeunits[2], '-', /EXTRACT)
        year=long(timeunits[0]);1850;long(STRMID(str[6],0,4))
        month=long(timeunits[1]);long(STRMID(str[6],4,2))
        day=long(timeunits[2])
    ;    if str[1] eq 'SIday' then begin
    ;      day=long(STRMID(str[6],6,2))
    ;    endif
        
        NCDF_CLOSE, file_ID
        print,'  All parameters of file ' + thesefiles[fidx] + ' have been writen out!'
        index=where(latitude[*,*] ge 47,count) ;只要纬度>47的数据
        ind = ARRAY_INDICES(latitude,index)
        for x=0,nt-1 do begin
          if timeinfo eq '360_day' then begin
            tempresult=calendar_360(year,month,day,time[x])
            yr=tempresult[0]
            mon=tempresult[1]
            dy=tempresult[2]
          endif else if timeinfo eq 'noleap' then begin
            tempresult=calendar_noleap(year,month,day,time[x])
            yr=tempresult[0]
            mon=tempresult[1]
            dy=tempresult[2]
          endif else begin
            caldat,julday(month,day,year)+time[x]-0.5,mon,dy,yr;thick需要在filename后面添加这个
    ;        print,mon,dy,yr
    ;        gregorian and julday For dates between 1 Jan CE and 15 Oct 1582 the two calendar systems differ by up to 10 days. For dates on or after 15 Oct 1582 the two calendar systems are identical.
          endelse
          
          
          if yr lt 1980 then begin
            continue
          endif
          outfilename = outfilepath +String(yr,FORMAT='(I04)')+'\'+strjoin(str[0:5],'_')+'_'+String(yr,FORMAT='(I04)')+String(mon,FORMAT='(I02)');+String(dy,FORMAT='(I02)');STRCOMPRESS(outfilepath +  'smos_sea_ice_thickness_'+year+'_'+month+'_'+day+'.tif',/remove_all)
          ;存储数据
          if str[1] eq 'SIday' then begin
            outfilename = outfilename+String(dy,FORMAT='(I02)')
          endif
          temp=[]
          if size(var_Data,/N_DIMENSIONS) eq 3 then begin ;这是由于变量不是像AWI，而是两维+1维时间
            temp=var_Data[*,*,x]
            temp=temp[index];var_Data[ind[1],ind[0],x]
          endif else begin
            temp=var_Data[index,x]
          endelse
          miss=where(temp eq varinfo,count)
          if count gt 0 then begin
            temp[miss]=-1
          endif
    ;      temp1=var_Data[ind[0],ind[1],x]
    ;      print,max(temp),mean(temp),var_Data[1234,137,x];,temp[100,37]
          dim=size(temp,/DIMENSIONS)
          if N_elements(dim) eq 1 then begin
            dim=[dim,1]
          endif
          envi_setup_head,FNAME=outfilename,NS=dim[0] ,NL=dim[1],NB=1,DATA_TYPE=size(temp,/type),INTERLEAVE=0,/OPEN,/WRITE ;4为float
          openw,lunou,outfilename,/get_lun
          writeu,lunou,temp;
          close,lunou
          free_lun,lunou
          ;消除内存
          temp=!Null
          delvar,temp
          print,outfilename
        endfor
        ;    输出经纬度
        dim=size(latitude[index],/DIMENSIONS)
        if N_elements(dim) eq 1 then begin
          dim=[dim,1]
        endif
        envi_setup_head,FNAME=infilepath+'\lat',NS=dim[0] ,NL=dim[1],NB=1,DATA_TYPE=size(latitude,/type),INTERLEAVE=0,/OPEN,/WRITE ;5为double
        openw,lunou,infilepath+'\lat',/get_lun
        writeu,lunou,latitude[index]
        close,lunou
        free_lun,lunou
        envi_setup_head,FNAME=infilepath+'\lon',NS=dim[0] ,NL=dim[1],NB=1,DATA_TYPE=size(longitude,/type),INTERLEAVE=0,/OPEN,/WRITE
        openw,lunou,infilepath+'\lon',/get_lun
        writeu,lunou,longitude[index]
        close,lunou
        free_lun,lunou
        ;    var_Data[where(~finite(var_Data))] = -1
        var_Data=!Null
        delvar,var_Data
        print,thesefiles[fidx],' is finished'
      ENDFOR
    
    
      nv.close
    
    
    end
    function calendar_360,yr,mon,dy,time
      tdy=dy+time mod 30
      tmon=long(time/30)
      tyr=yr+long(tmon/12)
      tmon=mon+tmon mod 12
      return,[tyr,tmon,tdy]
    end
    FUNCTION calendar_noleap,yr,mon,dy,time
      ;不考虑闰年
      ;
      months_days = [31,28,31,30,31,30,31,31,30,31,30,31]
      time=time+0.5
      ;  IF ISLEAPYEAR(year) EQ 1 THEN months_days[1]=28;29 ;注意我这边由于闰年的28都被排除了，所以DOY不计算29的情况
      tyr=yr+long(time/365)
      tdy=time mod 365
      for month=0,12-1 do begin
        if tdy le total(months_days[0:month]) then begin
          break
        endif
      endfor
      tmon=month+1
      if month gt 0 then begin
        tdy=tdy-total(months_days[0:month-1])
      endif
      
      RETURN, [tyr,tmon,tdy]
    END