up [pdf]
from rsf.proj import *

# download well-log data
Fetch('L-30.las','1406_Make_a_synthetic',
      server='https://raw.githubusercontent.com',
      top='seg/tutorials/master')

# Convert to RSF
Flow('L-30','L-30.las','las2rsf $SOURCE $TARGET',stdin=0,stdout=-1)

# Examine with "< L-30.rsf sfheaderattr segy=n desc=y"

# Mask bad values
Flow('mask1','L-30','headermath segy=n output=RHOB | mask min=0')
Flow('mask2','L-30','headermath segy=n output=DT   | mask min=0')
Flow('mask','mask1 mask2','mul ${SOURCES[1]}')

scale = dict(RHOB=1000,DT=3.28084,DEPTH=1/3.28084) # convert to SI units
threshold = dict(RHOB=100,DT=10)

for case in ('DT','RHOB','DEPTH'):
     # Extract values
     Flow(case,'L-30 mask',
          '''
          headermath segy=n output=%s | 
          headerwindow mask=${SOURCES[1]} | window |
          scale dscale=%g | 
          put o1=3058.5 d1=0.5 label1=Depth unit1=ft
          ''' % (case,scale[case]))

     # Clip spikes
     if case != 'DEPTH':
          despike = case+'-despike'
          Flow(despike,case,'despike wide=13')

          spikemask = case+'-spikemask'
          Flow(spikemask,[despike,case],
               '''
               math orig=${SOURCES[1]} output="abs(input-orig)" |
               mask min=%g
               ''' % threshold[case])

          Flow(despike+'2',[spikemask,despike,case],
               '''
               dd type=float |
               math orig=${SOURCES[2]} filt=${SOURCES[1]}
               output="input*filt+(1-input)*orig"
               ''')
          Result(case,[case,despike+'2'],
                 '''
                 cat axis=2 ${SOURCES[1]} |
                 graph title=%s
                 ''' % case)

# Compute time-to-depth relationship

log_start_time = 0.410553230106 # (s)

# windowed samples - 3837

Flow('dt1','L-30',
     '''
     headermath segy=n output=DT | 
     window n2=3837 | scale dscale=%g 
     ''' % scale['DT'])
Flow('dt2','dt1','window f1=21')
Flow('dt0','dt1 dt2','window n1=21 | math output=507.2 | cat axis=1 ${SOURCES[1]}')
Flow('deltat','dt0','stack axis=1 norm=n | scale dscale=%g' % (2*0.1524/1e6))

delta_time = 0.5596 # (s)

Flow('tdr','DT-despike','causint | math output="%g+(%g)*input" ' % (log_start_time+delta_time,2*0.1524/1e6))

Result('tdr','graph title="Twoway Time" yreverse=y transp=y label2=Time unit2=s')

# Compute acoustic impedance

Flow('z','RHOB slowness','math rho=${SOURCES[0]} slow=${SOURCES[1]} output="1e6*rho/slow" ')

# Compute reflection coefficient series

Flow('rc','z','ai2refl')

# Plots

Flow('dz','DEPTH','igrad')
Flow('slowness','DT-despike dz','div ${SOURCES[1]} | clip2 lower=300')

Plot('slowness','DEPTH slowness',
     '''
     cmplx ${SOURCES[1]} | window n1=21692 |
     graph transp=y yreverse=y min1=750 max1=4250 grid=y plotfat=3
     title=function unit1=m label2="P-wave slowness" unit2="\F10 m\F3 m/s" labelsz=12 titlesz=15
     ''')

Plot('twt','DEPTH tdr',
     '''
     cmplx ${SOURCES[1]} | window n1=21692 |
     graph transp=y yreverse=y min1=750 max1=4250 grid=y plotfat=3
     title=integral unit1=m label2="two-way time" unit2="s" labelsz=12 titlesz=15
     ''')

Plot('impedance','DEPTH z',
     '''
     cmplx ${SOURCES[1]} | window n1=21692 |
     graph transp=y yreverse=y min1=750 max1=4250 grid=y plotfat=3
     title=impedance unit1=m label2= unit2="kg/(m\^2\_s)" labelsz=12 titlesz=15
     ''')

Plot('reflectivity','DEPTH rc',
     '''
     cmplx ${SOURCES[1]} | window n1=21692 |
     graph transp=y yreverse=y min1=750 max1=4250 grid=y plotfat=3
     title=reflectivity unit1=m label2= unit2= labelsz=12 titlesz=15
     ''')

Plot('four','slowness twt impedance reflectivity','SideBySideAniso')

# add labels for well tops

welltops = {
    'Abenaki': 3404.311,
    'Base O-Marker': 2469.207,
    'Dawson Canyon': 984.504,
    'L Baccaro': 3964.534,
    'L Missisauga': 3190.646,
    'Logan Canyon': 1136.904,
    'Mid Baccaro': 3485.083,
    'U Missisauga': 2251.253,
    'Wyandot': 867.156
    }

a = (8.95-1.3)/(750-4250)
b = 1.3-4250*a
plots=[]
for top in welltops.keys():
    name = top.replace(' ','-').lower()
    Plot(name,None,'box font=2 x0=12.5 y0=%g label="%s" size=0.15 xt=0 yt=0' % (a*welltops[top]+b,top))
    plots.append(name)
Plot('welltops',plots,'Overlay')

Result('four','four welltops','Overlay')

Flow('tops.txt',None,
     'echo %s in=$TARGET n1=%d data_format=ascii_float' % (' '.join(map(str,welltops.values())),len(welltops)))
Flow('tops','tops.txt','dd form=native | sort')
Plot('tops','spray axis=1 n=2 o=0 d=1 | graph min2=750 max2=4250 pad=n yreverse=y wantaxis=n wanttitle=n')

Flow('ttops','tdr DEPTH tops',
     'iwarp warp=${SOURCES[1]} o1=932 d1=1 n1=3307 | inttest1 nw=4 interp=spline coord=${SOURCES[2]}')
Plot('ttops','spray axis=1 n=2 o=0 d=1 | graph min2=0.9 max2=2.9 pad=n yreverse=y wantaxis=n wanttitle=n')

# Convert to two-way-traveltime

Flow('zt','z tdr','cut f1=-20 | iwarp warp=${SOURCES[1]} o1=0 d1=0.004 n1=750')
Flow('rt','rc tdr','cut f1=-20 | iwarp warp=${SOURCES[1]} o1=0 d1=0.004 n1=750')

Plot('z','DEPTH z',
     '''
     cmplx ${SOURCES[1]} | window n1=21692 |
     graph transp=y yreverse=y min1=750 max1=4250 plotfat=3 plotcol=7
     title=impedance unit1=m label2= unit2="kg/(m\^2\_s)" labelsz=12 titlesz=15
     ''')
Plot('zw','z tops','Overlay')

Plot('zt',
     '''
     window min1=0.9 max1=2.9 |
     graph transp=y yreverse=y plotfat=3 plotcol=7
     title=impedance label1=Time unit1=s pad=n 
     label2= unit2="kg/(m\^2\_s)" labelsz=12 titlesz=15
     ''')
Plot('ztw','zt ttops','Overlay')

Plot('rt',
     '''
     window min1=0.9 max1=2.9 |
     graph transp=y yreverse=y plotfat=3 plotcol=7
     title=reflectivity label1=Time unit1=s pad=n 
     label2= unit2= labelsz=12 titlesz=15
     ''')
Plot('rtw','rt ttops','Overlay')

# Create synthetic

Flow('synth','rt','ricker1 frequency=25')
Plot('synth',
     '''
     window min1=0.9 max1=2.9 | 
     wiggle poly=y yreverse=y title=synthetic transp=y yreverse=y
     label1=Time unit1=s labelsz=12 titlesz=15 clip=0.002 grid=n plotcol=7
     ''')
Plot('synthw','synth ttops','Overlay')

Result('synth','zw ztw rtw synthw','SideBySideAniso')

# Get seismic data

Fetch('PenobXL_1155.txt','1406_Make_a_synthetic',
      server='https://raw.githubusercontent.com',
      top='seg/tutorials/master')

Flow('seismic','PenobXL_1155.txt',
     '''
     echo in=$SOURCE data_format=ascii_float n1=1003 n2=401 
     label2=Trace label1=Time unit1=s o1=0 d1=0.004 |
     dd form=native | reverse which=2 | window n1=750
     ''',stdin=0)

Plot('seismic','grey title=Seismic')
Flow('xpos',None,'spike n1=1 mag=77')
Plot('trace','synth xpos',
     '''
     wiggle transp=y yreverse=y wanttitle=n wantaxis=n plotcol=4
     xpos=${SOURCES[1]} xmin=400 xmax=0 poly=y zplot=10
     ''')
Result('seismic','seismic trace','Overlay')

# improve registration

End()

sflas2rsf
sfheadermath
sfmask
sfmul
sfheaderwindow
sfwindow
sfscale
sfput
sfdespike
sfmath
sfdd
sfcat
sfgraph
sfstack
sfcausint
sfai2refl
sfigrad
sfdiv
sfclip2
sfcmplx
sfbox
sfsort
sfspray
sfiwarp
sfinttest1
sfcut
sfricker1
sfwiggle
sfreverse
sfgrey
sfspike

https://raw.githubusercontent.com/seg/tutorials/master/1406_Make_a_synthetic/L-30.las
https://raw.githubusercontent.com/seg/tutorials/master/1406_Make_a_synthetic/PenobXL_1155.txt