一、获取CDF变量信息

pro test
  file = 'C:\Users\XCL\Desktop\omni_hro_1min_20080801_v01.cdf'
  ;打开cdf文件
  id = cdf_open(file)
  ;获取CDF文件全局变量信息
  inq = cdf_inquire(id)
  cdf_control,id,get_numattrs=num_attrs
  if (inq.nzvars ge 1) then begin
    for i=0, inq.nzvars-1 do begin
      var_inf=cdf_varinq(id,i,/zvariable)
      cdf_control,id,var=var_inf.name,get_var_info=info, /zvar
      if(info.maxrec le 0) then continue
      print,'-------------------------------'
      print,' ',var_inf.name
      print,'-------------------------------'
      for j=0, num_attrs[1]-1 do begin
        cdf_attinq,id,j+num_attrs[0],name1,scope1,maxentry1,maxzentry1
        if (cdf_attexists(id,j+num_attrs[0],i,/zvariable)) then begin
          cdf_attget,id,j+num_attrs[0],i,value1,/zvariable
          print,name1,value1,format='(2x,A10,":",4x,A)'
        endif
      endfor
    endfor
  endif
  print, 'Run over!'
end

得到基本参量信息后根据自己所需要的变量编写读取函数。本文需要的是磁场及太阳风相关

二、函数

pro read_cdf_OMNI,file,bx,by,bz,vx,vy,vz,pr,al
  id = cdf_open(file)
  a = cdf_inquire(id)
  cdf_control,id, variable='Epoch', get_var_info=info1
  cdf_control,id, variable='BY_GSM', get_var_info=info2
  cdf_control,id, variable='BZ_GSM', get_var_info=info3
  cdf_control,id, variable='BX_GSE', get_var_info=info4
  cdf_control,id, variable='Vx', get_var_info=info5
  cdf_control,id, variable='Vy', get_var_info=info6
  cdf_control,id, variable='Vz', get_var_info=info7
  cdf_control,id, variable='Pressure', get_var_info=info8
  cdf_control,id, variable='AL_INDEX', get_var_info=info9
  
  cdf_varget,id, 'Epoch', tt, rec_count=info1.maxrec, /zvariable
  cdf_varget,id, 'BY_GSM', by2, rec_count=info2.maxrec, /zvariable
  cdf_varget,id, 'BZ_GSM', bz2, rec_count=info3.maxrec, /zvariable
  cdf_varget,id, 'BX_GSE', bx2, rec_count=info4.maxrec, /zvariable
  cdf_varget,id, 'Vx', vx2, rec_count=info5.maxrec, /zvariable
  cdf_varget,id, 'Vy', vy2, rec_count=info6.maxrec, /zvariable
  cdf_varget,id, 'Vz', vz2, rec_count=info7.maxrec, /zvariable
  cdf_varget,id, 'Pressure', p2, rec_count=info8.maxrec, /zvariable
  cdf_varget,id, 'AL_INDEX', al2, rec_count=info9.maxrec, /zvariable
  cdf_close,id
  bx= {x: transpose(tt), y: transpose(bx2)}
  bad = where(bx.y gt 10,c)
  if c ne 0 then bx.y[bad] = !values.f_nan
  by = {x: transpose(tt), y: transpose(by2)}
  bad = where(by.y gt 10,c)
  if c ne 0 then by.y[bad] = !values.f_nan

  bz = {x: transpose(tt), y: transpose(bz2)}
  bad = where(bz.y gt 10,c)
  if c ne 0 then bz.y[bad] = !values.f_nan

  vx = {x: transpose(tt), y: transpose(vx2)}
  bad = where(vx.y gt 100,c)
  if c ne 0 then vx.y[bad] = !values.f_nan
  vy = {x: transpose(tt), y: transpose(vy2)}
  bad = where(vy.y gt 100,c)
  if c ne 0 then vy.y[bad] = !values.f_nan
  vz = {x: transpose(tt), y: transpose(vz2)}
  bad = where(vz.y gt 100,c)
  if c ne 0 then vz.y[bad] = !values.f_nan

  pr = {x: transpose(tt), y: transpose(p2)}
  bad = where(pr.y gt 2,c)
  if c ne 0 then pr.y[bad] = !values.f_nan
  al = {x: transpose(tt), y: transpose(al2)}
  bad = where(al.y gt 10,c)
  if c ne 0 then al.y[bad] = !values.f_nan
  return
end

三、绘图

pro last_hw
  read_cdf_OMNI,'C:\Users\XCL\Desktop\omni_hro_1min_20080801_v01.cdf',bx,by,bz,vx,vy,vz,pr,al
  t_om=CDF_EPOCH_TOJULDAYS(bx.x)
  
  trag=julday(8,25,2008,[6,10],0,0)
  ;设置图片尺寸相关参数
  x0=0.2;画图区域的左下角X坐标(归一化坐标系)
  xs=0.7;每一栏的长度
  y0=0.95;画图区域最大y坐标值
  ys=0.12;每一栏的高度
  yg=0.003;相邻两栏的间隔

  pngon,xsize=1800,ysize=2100
  loadct, 39
  plots, [0.491, 0.491], [0.34, 0.95], /normal, linestyle = 5, thick = 3,color=200
  plots, [0.551, 0.551], [0.34, 0.95], /normal, linestyle = 5, thick = 3,color=100
  plots, [0.2, 0.9], [0.89, 0.89], /normal, linestyle = 5, thick = 3
  plot,t_om,bx.y,/noerase,charsize=4,thick=4, $
    position=[x0,y0-1*ys-0*yg,x0+ xs,y0-0*ys-0*yg], $
    xthick=2,xticklen=0.1,xstyle=1,xminor=6,xtitle=' ',xtickinterval=3600/(24*60.0*60D),xtickname=replicate(' ',10), $
    ythick=2,yticklen=0.015,ystyle=1,xrange=trag,yrange=[-4,4],ytickinterval=2,yminor=4,ytitle=''
  oplot,t_om,by.y,thick=4,color=250;red 
  oplot,t_om,bz.y,thick=4,color=50
  axis,x0-0.04, y0 - 1*ys - 0*yg, yaxis=0, yticklayout=1, /normal, ytickname =replicate(' ',10), charsize=5,ytitle = 'B!dIMF!n!c[nT]'
  
  plot,t_om,vx.y,/noerase,charsize=4,thick=4, $
    position=[x0,y0-2*ys-1*yg,x0+ xs,y0-1*ys-1*yg], $
    xthick=2,xticklen=0.1,xstyle=1,xminor=6,xtitle=' ',xtickinterval=3600/(24*60.0*60D),xtickname=replicate(' ',10), $
    ythick=2,yticklen=0.015,ystyle=1,xrange=trag,yrange=[-304,-280],ytickinterval=10,yminor=5,ytitle=''
  axis,x0-0.04, y0 - 2*ys - 1*yg, yaxis=0, yticklayout=1, /normal, ytickname =replicate(' ',10),charsize=5, ytitle = 'V!dX,SW!n!c[km/s]'
  
  plot,t_om,vy.y,/noerase,charsize=4,thick=4, $
    position=[x0,y0-3*ys-2*yg,x0+ xs,y0-2*ys-2*yg], $
    xthick=2,xticklen=0.1,xstyle=1,xminor=6,xtitle=' ',xtickinterval=3600/(24*60.0*60D),xtickname=replicate(' ',10), $
    ythick=2,yticklen=0.015,ystyle=1,xrange=trag,yrange=[-20,4],ytickinterval=10,yminor=5,ytitle=''
  oplot,t_om,vz.y,thick=4,color=250
  axis,x0-0.04, y0 - 3*ys - 2*yg, yaxis=0, yticklayout=1, /normal, ytickname =replicate(' ',10),charsize=5, ytitle = 'V!dYZ,SW!n!c[km/s]'
 
  plot,t_om,pr.y,/noerase,charsize=4,thick=4, $
    position=[x0,y0-4*ys-3*yg,x0+ xs,y0-3*ys-3*yg], $
    xthick=2,xticklen=0.1,xstyle=1,xminor=6,xtitle=' ',xtickinterval=3600/(24*60.0*60D),xtickname=replicate(' ',10), $
    ythick=2,yticklen=0.015,ystyle=1,xrange=trag,yrange=[0.7,1.2],ytickinterval=0.2,yminor=2,ytitle=''
  axis,x0-0.04, y0 - 4*ys - 3*yg, yaxis=0, yticklayout=1, /normal, ytickname =replicate(' ',10), charsize=5,ytitle = 'P!dSW!n!c[nPa]'
  
  plot,t_om,al.y,/noerase,charsize=4,thick=4, $
    position=[x0,y0-5*ys-4*yg,x0+ xs,y0-4*ys-4*yg], $
    xthick=2,xticklen=0.1,xstyle=1,xminor=6,xtitle=' ',xtickinterval=3600/(24*60.0*60D),xtickname=replicate(' ',10), $
    ythick=2,yticklen=0.015,ystyle=1,xrange=trag,yrange=[-150,0],ytickinterval=50,yminor=5,ytitle=''
  axis,x0-0.04, y0 - 5*ys - 4*yg, yaxis=0, yticklayout=1, /normal, ytickname =replicate(' ',10), charsize=5,ytitle = 'AL!n!c[nT]' 

  xyouts,0.74,0.96,'2008-08-25',charsize = 6,CHARTHICK = 1,/normal
  xyouts,0.2,0.96,'OMNI',charsize = 6,CHARTHICK = 1,/normal
  xyouts, 0.91, 0.93, 'B!dx', charsize = 5, /normal, color = 0
  xyouts, 0.91, 0.905, 'B!dy', charsize = 5, /normal, color = 250
  xyouts, 0.91, 0.88, 'B!dz', charsize = 5, /normal, color = 50  
  xyouts,0.1,0.31,'UT',charsize = 5,charthick =5,/normal
  xyouts,0.17,0.31,'0600',charsize = 5,/normal
  xyouts,0.35,0.31,'0700',charsize = 5,/normal
  xyouts,0.53,0.31,'0800',charsize = 5,/normal
  xyouts,0.71,0.31,'0900',charsize = 5,/normal
  xyouts,0.88,0.31,'1000',charsize = 5,/normal
 
  xyouts, 0.91, 0.68,'V!dy',charsize = 5, /normal, color = 0
  xyouts, 0.91, 0.655,'V!dz',charsize = 5, /normal, color = 250

  pngoff,'E:\Plot3.png'

end

四、运行效果

补充: 此代码需要安装spedas库。pngon、pngoff两个函数上篇文章已记录

 

更多推荐