IDL 绘制线形图

2013-10-29 16:20

今天尝试用IDL做了一张图，以前在IDL中也做过类似的工作，但今天是用的Using IDL in Meteorology的工具包，并且保存为ps格式的文件。

先上图：

 其实做来做去，基本呈现结果都一样。但通过这样的工作，在利用直接图形法的应用上，感觉不一样了。以前觉得快速可视化是基于对象图形法的，用起来非常方便，现在看来，其实直接图形法也很强大，而且直接图形法是IDL制图的最基础的内容。

总的来看，有以下几点：

1）想了个办法去画这个阴影。以前是直接用POLYGON命令来实现的，现在用直接图形法是怎么实现的呢？其实是转了个弯来实现这个阴影绘制的。就是先用上下两个边界上所有的点，来构成一个多边形（而且是闭合的），然后用一个gplot中的fillcol参数即可实现阴影。

2）直接图形法做图例还是很欠缺的，没办法进行复杂的图例设置。
==============================CODE===========================

pro test
  CD, 'C:\Users\WangKang\IDLWorkspace80\Default'
 
  ;   Read into the data file ===
 
  data =READ_ASCII('change_Min.nums.dat')
  data =data.field1
 
  year = data(0,*)
 
  avr = data(1,*)
 
  std = data(2,*)
 
  y_up = avr+std ;Up boundary
 
  y_dn = avr-std ;Low boundary
 
  ; Draw figure and save in ps file type ===
 
  PSOPEN $
    ,AXISTYPE=2 $
    ,XSIZE=21000 $
    ,YSIZE=12000 $
    ,YOFFSET=2500 $
    ,file='idl.ps'
   
  CS,SCALE=7,NCOLS=13
 
  GSET,xmin=min(year),xmax=max(year),ymin=-30,ymax=20
 
  ; Make up a poly
 
  x0= [[year], [REVERSE(year,2)], [year[0]]]
  y0= [[y_up], [REVERSE(y_dn,2)], [y_up[0]]]

x0= REFORM(x0,n_elements(x0))
  y0= REFORM(y0,n_elements(y0))
 
  ; Fill the ploy
 
  gplot,x=x0,y=y0,FILLCOL=5,/NOLINES
 
  ; Plot some lines on the filled poly
 
  gplot $
    ,col=[1,1] $
    ,STYLE=[5, 5] $
    ,LEGPOS=3 $
    ;    ,/LEGEND $
    ,LABELS=['Up','Down'] $
    ,x=transpose(year),y=transpose([y_up,y_dn])
   
  ; Low pass filter
   
  cutoff=0.11
 
  order = 1
 
  y_low_pass = low_pass(avr,cutoff,order)
 
  CS,SCALE=2,NCOLS=22
 
  gplot $
    ,col=[1,16] $
    ,STYLE=[0, 0] $
    ,LEGPOS=3 $
    ,/LEGEND $
    ,/NOLEGENDBOX $
    ,LABELS=['Mean','Low Pass Filter'] $
    ,x=transpose(year),y=([[avr],[y_low_pass]])

   
  ; Regression estimate ===
   
  coef=imsl_multiregress(transpose(year),avr,anova_table=anova);

p_value=anova[9];
 
  gplot,x=year,y=year*coef(1)+coef(0), STYLE=2
 
  GPLOT, x=1984, y =-28, text='Cut off='+string(cutoff,FORMAT='(F5.2)')+ $
 
    '; Slope='+string(coef(1),FORMAT='(F5.2)') + '; P='+string(p_value,FORMAT='(F5.3)')
   
  ; Make up the axes
   
  axes, XTITLE='Year',YTITLE='Number of Frozen Days',XMINOR=1,YMINOR=1
 
  PSCLOSE
 
end

function low_pass,s,cutoff,order
  s0=s;
 
  s = s0(where(finite(s0)));
  n=size(s,dimension=1);
 
  n=n[0];
 
  s_hat = fft(s);
 
  t=indgen(n)+1;
  dt=t[1]-t[0];
 
  c = linfit(t, s)
 
  s_fit = c[0] + c[1]*t
 
  s_dtr = s - s_fit

  welch = 1 - (findgen(n) / (n-1)*2 - 1)^2
 
  s_win = s_dtr * welch
 
  freq_c = cutoff
 
  kernel = (dist(n))[*,0]/(n*dt)
  filter = 1 /(1 + (kernel/freq_c)^(2*order))
  s_hat = fft(s_dtr)
  s_f = fft(s_hat*filter, /inverse)
  s_f = s_f + c[0] + c[1]*t;
 
  return,s_f
 
end

 

===============================================================

返回