IDL CMIP6 NC格式数据处理

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

AWI-CM-1-1-MR
CMCC.CMCC-CM2-HR4
CNRM-CERFACS.CNRM-CM6-1-HR
NOAA-GFDL.GFDL-CM4
MOHC.HadGEM3-GC31-MM

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

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

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

最后也保存了lon和lat

pro read_CMIP6_model_data
  infilepath = 'H:\Model\AWI-CM-1-1-MR\historical\scon\daily'
  outfilepath = 'H:\Model\AWI-CM-1-1-MR\historical\scon\daily\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('*.nc',count = numFiles)
  print,'number of files found:',numFiles
  ;处理数据
  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)
    ;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 latid eq -1 then begin
      latid = 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        ;获取变量数据
    ;大于等于0为海冰有效值,缺失值为-0.2,陆地为-0.1


    year=long(STRMID(str[6],0,4))
    month=long(STRMID(str[6],4,2))
    day=1
    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的数据
    for x=0,nt-1 do begin
      caldat,julday(month,day,year)+time[x]-0.5,mon,dy,yr;thick需要在filename后面添加这个
      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=var_Data[index,x]

      envi_setup_head,FNAME=outfilename,NS=count,NL=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
      ;消除内存

      delvar,temp
      print,outfilename
    endfor
    ;    输出经纬度
    envi_setup_head,FNAME=infilepath+'\lat',NS=count,NL=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=count,NL=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

    print,thesefiles[fidx],' is finished'
  ENDFOR


  nv.close


end

 

上一篇:微信小程序实现定位功能


下一篇:微信小程序简单封装获取定位