波形图

比如，绘制子波形态，单道地震记录等。采用sfgraph来绘图，类似与Matlab的plot命令。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

# plot wavelet
from rsf.proj import *

par={
}

def waveplot(custom,par):
    return '''
    graph title="" min1=%g min2=0 max2=+1
    plotfat=5 plotcol=5
    label1=%s unit1=%s
    label2="Amplitude" unit2=""
    parallel2=n screenratio=0.5 screenht=7
    %s
    ''' % (par['oxx'],par['lxx'],par['uxx'],
           par['labelattr']+' '+custom)

par['labelattr']=' parallel2=n labelsz=6 labelfat=3 titlesz=12 titlefat=3 '

par['lxx']='x'
par['uxx']='m'
par['oxx']= 0.

Result('wav','transp |' + waveplot('',par))

说明wav为绘图的数据，其维度必须为1xn2的形式。

示例：暂缺

3D图

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

par = {
    'nx':41, 'ox':0, 'dx':0.05,  'lx':'x', 'ux':'km',
    'ny':41, 'oy':0, 'dy':0.05,  'ly':'y', 'uy':'km',
    'nz':41, 'oz':0, 'dz':0.05,  'lz':'z', 'uz':'km'
    }
# 生成三维数据
Flow('data',None,'''math n1=41 n2=41 n3=41 o1=-1 o2=-1 o3=-1
     d1=0.05 d2=0.05 d3=0.05 output="exp(-2*(x1^2+x2^2+x3^2))"''')
# 绘图
fdmod.param(par)
# 立体图
Result('data', 'byte gainpanel=a pclip=95 |'
+fdmod.cgrey3d('flat=y ',par))
# 三视图
Result('data1','data','byte gainpanel=a pclip=95 |'
+fdmod.cgrey3d('flat=n ',par))

立体图

立体图.png

三视图

三视图.png

绘制box

常用于论文中的局部比较

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

from rsf.proj import *
from rsf.recipes import fdmod,rwezo
from math import pi

zmin=0
zmax=1
xmin=0
xmax=1

# make box data
fdmod.makebox('box',zmin,zmax,xmin,xmax,par)

#           __________________                      zmin
#          |                  |
#          |                  |
#          |                  |
#          |                  |
# xmin     |__________________|  xmax               zmax

# draw box
Plot('box',fdmod.bbplot('',par))

box.png

半圆

常用于脉冲响应的绘制。

给定参数

``` {.py}
  xcenter = 1285.0　# 圆心x坐标
  zcenter = 0.0　　 # 圆心z坐标　
  sampling = 1000　# 采样点数
  radius = 768　　　# 半径
```

绘图

``` {.py}
  fdmod.circle('cc',xcenter,zcenter,radius,sampling,par)
  Plot(  'cc',fdmod.ssplot('plotfat=7 plotcol=2 ',par))
```

结果

半圆图.png

椭圆

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

  #定义临时处理Flow
  def Temp(o,i,r):
    Flow(o,i,r+ ' datapath=%s '%os.environ.get('TMPDATAPATH',os.environ.get('DATAPATH')))
　#椭圆模块
  def ellipse(cc,xcenter,zcenter,semiA,semiB,sampling,par):
    Temp(cc+'_x',None,
         'math n1=%d d1=%g o1=%g output="%g+%g*cos(%g*x1/180)"'
         % (sampling,360./sampling,0.,xcenter,semiA,pi) )
    Temp(cc+'_z',None,
         'math n1=%d d1=%g o1=%g output="%g+%g*sin(%g*x1/180)"'
         % (sampling,360./sampling,0.,zcenter,semiB,pi) )
    Flow(cc,[cc+'_x',cc+'_z'],
         '''
         cat axis=2 space=n
         ${SOURCES[0]} ${SOURCES[1]} | transp |
         put label1="" unit1="" label2="" unit2=""
         ''', stdin=0)

xcenter=(par['ox']+(par['nx']-1)*par['dx'])/2.0　# 坐标系中心x
zcenter=par['oz']  ＃坐标系中心z
ab=par['dx']*(361-150)# 计算焦距
semiA=(par['kt']-50)*par['dt']*1500　# 计算a
semiB=np.sqrt(semiA**2.0-(ab/2.0)**2.0) # 计算b
ellipse('ec',xcenter,zcenter,semiA,semiB,500,par) # 生成椭圆坐标

Plot(  'ec',fdmod.ssplot('plotfat=10 plotcol=2 symbol=.',par))　# 绘制椭圆

结果

椭圆.png

椭圆坐标系

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

# 生成坐标系文件 rays.rsf
Flow('rays','datRC velo',
   '''
   %(rays)slrays
	vel=${SOURCES[1]}
	dt=%(dt)g nt=%(nt)d
	ot=%(ot)g minang=%(minang)g maxang=%(maxang)g
 '''%par)
# 绘制椭圆坐标系
rwezo.cos('rays',20,100,'',par)

椭圆坐标系.png