from rsf.proj import *

graph = 'graph plotcol=1 wanttitle=n min2=-6 max2=6 min1=-1 max1=19 screenwd=12 screenht=4 labelsz=6 '

# Make a trace
Flow('trace',None,'echo 0,1,3,0,-2,-5,-2,2,5,0 | csv2rsf | put o1=0 d1=2 label1= unit1=')
Plot('trace',graph + 'symbol=o symbolsz=5 plotcol=1')

# Linear interpolation
Flow('linear','trace','remap1 n1=901 d1=0.02 order=2')
Plot('linear',graph + 'plotcol=1')

Flow('pi',None,'math n1=1 output=3.14159')
Flow('fpi','trace pi','inttest1 coord=${SOURCES[1]} interp=lag nw=2')

Flow('pi2','pi','spray axis=1 n=2')
Flow('pi0','pi','math output=-1 | cat axis=1 $SOURCE')

Plot('pi','fpi pi2',
     '''
     math output=-6 | cat axis=1 ${SOURCES[0]} | cat axis=2 ${SOURCES[1]} order=2,1 |
     transp | dd type=complex | window | %s plotcol=2
     ''' % graph)

Plot('fpi','fpi pi0',
     '''
     spray axis=1 n=2 | cat axis=2 ${SOURCES[1]} order=2,1 |
     transp | dd type=complex | window | %s plotcol=2
     ''' % graph)

Plot('pifpi','pi fpi','cmplx ${SOURCES[1]} | %s symbol=o symbolsz=5 plotcol=2' % graph)

Result('linear','linear trace pi fpi pifpi','Overlay')

graph += 'max1=20 '

# Make a trace
Flow('trace2',None,
     '''
     echo 0,1,3,0,-2,-5,-2,2,5,0,-1,1,0,3,3,-1,-3,-5,-1,0 | 
     csv2rsf | put n2=1 o2=0 o1=0 d1=1 label1= unit1=
     ''')
Plot('trace2',graph + 'symbol=o symbolsz=5 plotcol=1')

grid = 'n1=2000 o1=0 d1=%g' % (19.0/1999)

Flow('spline','trace2','spline %s' % grid)
Plot('spline',graph + ' plotcol=2')

Flow('linear2','trace2','remap1 %s order=2' % grid)
Plot('linear2',graph + ' plotcol=4')

Flow('x','trace2','math output=x1 | transp | pad n1=2')
Flow('nearest','trace2 x','nnint coord=${SOURCES[1]} %s voro=y n2=1 o2=0 d2=1' % grid)
Plot('nearest',graph + ' plotcol=6 dash=2')

Result('trace2','nearest linear2 spline trace2','Overlay')

# Real data

sgy = 'Penobscot_0-1000ms.sgy'

Fetch(sgy+'.gz','agilegeo',
      server='https://s3.amazonaws.com',top='')
zcat = WhereIs('gzcat') or WhereIs('zcat')
Flow(sgy,sgy+'.gz',zcat + ' $SOURCE',stdin=0)
Flow('seismic tseismic',sgy,'segyread tfile=${TARGETS[1]}')

Flow('cube','seismic','intbin xk=iline yk=xline')

Result('seismic','cube',
       '''
       window n2=1 f2=99 | 
       grey title=Seismic scalebar=y screenht=9 screenwd=18 barlabel=Amplitude
       ''')

# Horizon

Fetch('Penobscot_Seabed.txt','data',
      server='https://raw.githubusercontent.com',
      top='seg/tutorials-2016/master/1604_Function_of_interpolation')
Flow('xy','Penobscot_Seabed.txt',
     '''
     echo in=$SOURCE data_format=ascii_float n1=3 n2=249572 | 
     dd form=native type=int | window n1=2
     ''',stdin=0)
Flow('horizon','Penobscot_Seabed.txt xy',
     '''
     echo in=$SOURCE data_format=ascii_float n1=3 n2=249572 | dd form=native | 
     window n1=1 f1=2 squeeze=n | intbin head=${SOURCES[1]} xkey=0 ykey=1 xmin=1000 xmax=1600 ymin=1000 ymax=1480 nx=601 ny=481 |
     window |  scale dscale=0.001 | put label1=Inline label2=Crossline
     ''',stdin=0)

Result('horizon','grey color=gist_earth scalebar=y title=Horizon yreverse=n allpos=y bias=0.015 polarity=y')

# Extract amplitude on one line

line=300-1

Flow('horizon2','horizon','window n1=1 f1=%d' % line)
Plot('horizon2',
     '''
     window min1=1003 max1=1477 |
     graph min2=0 max2=1 min1=1000 max1=1480 yreverse=y scalebar=y screenht=9 
     screenwd=18 plotcol=3 plotfat=5 wanttitle=n wantaxis=n
     ''')

Plot('seismic2','cube',
     '''
     window n2=1 f2=%d | 
     grey wanttitle=n scalebar=y screenht=9 screenwd=18 barlabel=Amplitude
     ''' % line)

Result('seismic2','seismic2 horizon2','Overlay')

# Zoom in

Flow('horizon2-zoom','horizon2','window n1=150')
Plot('horizon2-zoom',
     '''
     window min1=1003 |
     graph min2=0.1 max2=0.3 min1=1000 max1=1151 yreverse=y scalebar=y screenht=9 
     screenwd=18 plotcol=3 plotfat=10 wanttitle=n wantaxis=n
     ''')

Flow('seismic2-zoom','cube','window n2=1 f2=%d n3=150 min1=0.1 max1=0.3' % line)
Plot('seismic2-zoom',
     '''
     grey wanttitle=n scalebar=y screenht=9 screenwd=18 barlabel=Amplitude
     ''')

Result('seismic2-zoom','seismic2-zoom horizon2-zoom','Overlay')

# Extract amplitude

Flow('coord-zoom','horizon2-zoom','transp')
Flow('ampl-zoom','seismic2-zoom coord-zoom',
     '''
     inttest1 coord=${SOURCES[1]} same=n interp=spline nw=4 | window
     ''')
Result('ampl-zoom',
       '''
       window min1=1003 | 
       graph title="Extracted Amplitude" unit1= unit2= label2=Amplitude
       ''')

# Amplitude of entire horizon

Flow('coord','horizon','spray axis=1 n=1')
Flow('ampl','cube coord',
     '''
     inttest1 coord=${SOURCES[1]} same=n interp=spline nw=4 | window
     ''')
Result('ampl','grey title=Amplitude unit1= yreverse=n polarity=y allpos=y')

# Compare different methods

methods = {'spl': 'interp=spline nw=4',
           'near': 'interp=lag nw=1',
           'lin': 'interp=lag nw=2'}
for m in methods.keys():
    Flow(m,'seismic2-zoom coord-zoom',
         '''
         inttest1 coord=${SOURCES[1]} same=n %s | window
         ''' % methods[m])
Result('methods','spl lin near',
       '''
       cat axis=2 ${SOURCES[1:3]} | window min1=1003 | 
       graph title="Extracted Amplitude" unit1= unit2= label2=Amplitude
       dash=0,1,2
       ''')

# Timeslice half-way through the data

Flow('half',None,'math n1=1 output=%g' % (124.5*0.004))
Flow('seismic2','cube','window n2=1 f2=%d n3=150' % line)

for m in methods.keys():
    Flow(m+'-half','seismic2 half',
         '''
         inttest1 coord=${SOURCES[1]} same=y %s | window
         ''' % methods[m])
Result('half','spl-half lin-half near-half',
       '''
       cat axis=2 ${SOURCES[1:3]} | window min1=1003 | 
       graph title="Extracted Amplitude" unit1= unit2= label2=Amplitude
       dash=0,1,2
       ''')

# Amplitudes along a deviated well-bore

Flow('trajectory.txt',None,
     '''echo   
     0   0    0
     0   0 -100
     0   0 -200
     5   0 -300
    10  10 -400
    20  20 -500
    40  80 -650
   160 160 -700
   600 400 -800
  1500 960 -800
     n1=3 n2=10 data_format=ascii_int 
     in=$TARGET
     ''')
Flow('trajectory','trajectory.txt',
     'dd form=native type=float | transp')

# Shift the 'well' to its tophole location.
Flow('txy','trajectory','window n2=2 | add add=3000')
Flow('trajectory2','trajectory txy',
     '''
     window f2=2 | 
     cat ${SOURCES[1]} axis=2 order=2,1 |
     put o1=0 d1=1 |
     spline n1=1000 o1=0 d1=0.009
     ''')

Flow('x2','trajectory2','window n2=1 f2=0')
Flow('y2','trajectory2','window n2=1 f2=1')
Flow('z2','trajectory2','window n2=1 f2=2')

Plot('xy','x2 y2','cmplx ${SOURCES[1]} | graph title=XY')
Plot('xz','x2 z2','cmplx ${SOURCES[1]} | graph title=XZ')
Plot('yz','y2 z2','cmplx ${SOURCES[1]} | graph title=YZ')
Result('trajectory','xy xz yz','SideBySideAniso')

Flow('cubet','cube','transp plane=13 memsize=1000 | put o1=0 d1=25 o2=0 d2=12.5 o3=0 d3=-4')
Flow('trajectoryt','trajectory2','transp')

Flow('trace3','cubet trajectoryt',
     'inttest3 coord=${SOURCES[1]} interp=lag nw=2 | put unit1= label1=Sample d1=1 o1=1')
Result('trace3','wiggle poly=y title="Deviated Well" pclip=100')

End()