Merge pull request #7 from Dirack/feature/issue/2

Feature/issue/2

experiments/README.md

@@ -0,0 +1 @@
+#### Experiments and studies of CRE Velocity Inversion

experiments/algorithmVelocityInversion/README.md

@@ -0,0 +1 @@
+#### Scratch of velocity inversion algorithm

experiments/algorithmVelocityInversion/SConstruct

@@ -0,0 +1,191 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# SConstruct  (Madagascar Script)
+#
+# Purpose: Scratch of velocity inversion algorithm.
+#
+# Site: http://www.dirackslounge.online
+# 
+# Version 1.0
+#
+# Programer: Rodolfo A. C. Neves (Dirack) 19/08/2020
+#
+# Email: rodolfo_profissional@hotmail.com
+#
+# License: GPL-3.0 <https://www.gnu.org/licenses/gpl-3.0.txt>.
+
+# Madagascar package
+from rsf.proj import *
+
+# Personal Madagascar recipes
+from rsf.recipes.kimodel import multiLayerModel as mlmod
+from rsf.recipes.kimodel import kirchhoffNewtonModeling as kinewmod
+from rsf.recipes.velocityAnalysis import velocityAnalysis as nmoStack
+from rsf.recipes.diff import diffmig as diff
+
+# Velocity continuation recipe
+from rsf.recipes.velcon import velcon
+
+xmax = 6.0
+zmax = 2.0
+
+layers = ((0.30,0.50,0.20,0.30),
+	  (1.65,1.85,1.55,1.65))
+
+velocities = (1.508,
+	      1.690,
+	      2.0)
+
+# Generate multi layer model and data cube
+mlmod(interfaces='interfaces',
+	dipsfile='interfacesDip',
+	modelfile='mod1',
+	xmax=xmax,
+	zmax=zmax,
+	layers=layers,
+	velocities=velocities)
+
+kinewmod(reflectors='interfaces',
+	reflectorsDip='interfacesDip',
+	filename='multiLayerDataCube',
+	velocities=velocities,
+	nt=1001,
+	dt=0.004,
+	ns=201,
+	ds=0.025,
+	nh=161,
+	dh=0.025)
+
+Flow('dataTransposed','multiLayerDataCube','transp plane=23 | transp plane=34')
+
+# Pre-stack velocity continuation
+velcon(data='dataTransposed',        # data name
+           nv=100,          # continuation steps
+           v0=1.5,          # initial velocity
+           dv=0.01,          # velocity step
+           nx=201,          # lateral dimension
+           nh=161,          # number of offsets
+           padt=1024,        # time padding
+           padt2=2048,       # extra time padding
+           padx=None,   # lateral padding
+           v1=None,     # other velocity
+           n1=None,     # time extent
+           dt=0.004,    # time sampling
+           dx=0.025,     # midpoint sampling
+           units='km',  # lateral units
+           vslope=None, # semblance muting
+           vx0=0,       # semblance muting
+           x0=0,        # lateral origin
+           srect1=3,    # semblance vertical smoothing
+           srect2=1,    # semblance lateral smoothing
+           rect1=10,    # vertical smoothing
+           rect2=10)   # lateral  smoothing
+
+# velocity analysis, NMO correction and stack
+nmoStack(dataCube='multiLayerDataCube',
+	pick='vnmo',
+	stack='stackedSection',
+	vrms='vrms',
+	v0=1.5,
+	dv=0.01,
+	nv=100,
+	vel0=1.5,
+	rect1=15,
+	rect2=40,
+	rect3=3,
+	dt=0.004)
+
+# Post-stack velocity continuation
+Flow('velocityCube','stackedSection',
+	'''
+	pad beg2=200 end2=200 | cosft sign2=1 |
+	stolt vel=1.5 |
+	vczo v0=1.5 dv=0.01 nv=100 verb=y |
+	transp plane=23 | cosft sign2=-1 |
+	window min2=0 max2=5
+	''')
+
+# Loop over reflectors
+numberOfReflectors = len(layers)
+section = 'stackedSection'
+
+for i in range(numberOfReflectors):
+
+	reflectorPickedPoints = 'reflectorPickedPoints-%i.txt' % i
+	t0sFile = 't0s-%i' % i
+	t0sAscii = 't0s-%i.asc' %i
+	m0sFile = 'm0s-%i' % i
+	m0sAscii = 'm0s-%i.asc' % i
+	returnedSection = 'returnedSection-%i' % i
+	diffSection = 'diffSection-%i' % i
+
+	# Reflector iterative Picking
+	Flow(reflectorPickedPoints,section,
+		'''
+		ipick
+		''')
+
+	# Next step is done with 'ascFormat.sh' Shell Script
+	# please check if the script have permissions to execute
+	# Build t0 coordinates file (pass 1 to generate t0s file)
+	Flow(t0sAscii,reflectorPickedPoints,
+		'''
+		./ascFormat.sh 1 %s
+		''' % (t0sAscii))
+
+	Flow(t0sFile,t0sAscii,'sfdd form=native')
+
+	# Build m0 coordinates file (pass 2 to generate m0s file)
+	Flow(m0sAscii,reflectorPickedPoints,
+		'''
+		./ascFormat.sh 2 %s
+		''' % (m0sAscii))
+
+	Flow(m0sFile,m0sAscii,'sfdd form=native')
+
+	# Diffraction simulation in stacked section
+	Flow([returnedSection,diffSection],[section,t0sFile,m0sFile],
+		'''
+		diffsim diff=${TARGETS[1]} aperture=1
+		t0=${SOURCES[1]} m0=${SOURCES[2]} v=%g freq=10 verb=y
+		''' % (velocities[i]))
+
+	section = returnedSection
+
+Flow('diffAdd',['diffSection-0','diffSection-1'],'sfadd ${SOURCES[1]} scale=1,1')
+
+diff("diff",
+	'diffAdd',
+	v0=1.4,
+	nv=100,
+	dv=0.01,
+	nx=201,
+	padx=1000,
+	nt=1000,
+	tmin=0,
+	tmax=4,
+            rect1=10,
+            rect2=10,
+            srect1=1,
+            srect2=3,
+            vslope=None,
+            units='Km',
+            f1=1,
+            j3=1,
+            dx=0.025,
+            x0=0,
+            beg1=0,
+            frect1=0,
+            frect2=0,
+            an=1,
+            nout=2048,
+            vx0=None)
+
+# Use aliases to split building
+# Call it 'scons alias' to build
+# target identified by alias
+Alias('model','multiLayerDataCube.rsf')
+Alias('stack','velocityCube.rsf')
+
+End()

experiments/algorithmVelocityInversion/ascFormat.sh

@@ -0,0 +1,35 @@
+#!/bin/bash
+#
+# ascFormat.sh (Shell Script)
+# 
+# Purpose: Convert sfipick picked points
+# to Madagascar readable ascii format.
+# Important! $1=1 cut and format the first collumn of the file (t0 time)
+# If $1=2 cut and format the second collumn of the file (m0 CMP)
+# 
+# Site: https://dirack.github.io
+# 
+# Version 1.0
+# 
+# Programer: Rodolfo A C Neves (Dirack) 20/08/2020
+# 
+# Email: rodolfo_profissional@hotmail.com
+# 
+# License: GPL-3.0 <https://www.gnu.org/licenses/gpl-3.0.txt>.
+
+INPUT=$(cat "-" | tr '\t' ' ' | cut -d' ' -f"$1")
+NUMBER_OF_POINTS=$(echo "$INPUT" | wc -l)
+FORMATED_INPUT=$(echo "$INPUT" | tr '\n' ' ')
+
+if [ "$1" == "1" ]
+then
+	FILE="t0s"
+elif [ "$1" == "2" ]
+then
+	FILE="m0s"
+else
+	echo "Error: $(basename $0): \$1 should be equal to 1 or 2!"
+	exit 1
+fi
+
+echo "${FORMATED_INPUT}n1=$NUMBER_OF_POINTS d1=1 o1=0 data_format=ascii_float in=$2"

experiments/modelUpdateStudy/README.md

@@ -0,0 +1 @@
+#### Study about model update and interface interpolation with cubic spline

experiments/modelUpdateStudy/SConstruct

@@ -0,0 +1,110 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# SConstruct  (Madagascar Script)
+#
+# Purpose: Study of interface interpolation with cubic spline
+# and velocity model update.
+#
+# Site: http://www.dirackslounge.online
+# 
+# Version 1.0
+#
+# Programer: Rodolfo A. C. Neves (Dirack) 07/06/2020
+#
+# Email: rodolfo_profissional@hotmail.com
+#
+# License: GPL-3.0 <https://www.gnu.org/licenses/gpl-3.0.txt>.
+
+# Madagascar package
+from rsf.proj import *
+
+# Error report function
+import atexit
+from errorReport import print_build_failures
+
+from random import uniform, seed
+
+seed()
+xmax = 6.0
+zmax = 2.0
+
+def arr2str(array,sep=' '):
+	'''
+	Convert a tuple into a comma separeted string
+	'''
+	return string.join(map(str,array),sep)
+
+# Velocity Model update
+for i in range(1):
+
+	interfaces = 'iterfaces-%i' %i
+	dipsfile = 'dipsfile-%i' %i
+	modelfile = 'model-%i' %i
+	ascInterfaces = 'layers-%i.asc' %i
+	layersFile = 'layers-%i' %i 
+	raysfile = 'rays-%i' %i
+	huygens = 'huygens-%i' %i
+	overlay = 'overlay-%i' %i
+
+	layers = [0.65,0.85,0.55,0.65]
+	layers = map(lambda x: x+uniform(-0.1,0.1),layers)
+
+	velocities = [1.508,2.0]
+
+	vstr=arr2str(velocities,',')
+
+	n1 = 4 #len(layers[0])
+	n2 = 1 #len(layers)
+
+	Flow(ascInterfaces,None,
+	     '''
+	     echo %s
+	     n1=%d n2=%d o1=0 d1=%g
+	     data_format=ascii_float in=$TARGET
+	     ''' % (arr2str(layers),n1,n2,xmax/(n1-1)))
+
+	Flow(layersFile,ascInterfaces,'dd form=native')
+
+	d = 0.00101 # non-round for reproducibility
+
+	Flow(interfaces,layersFile,
+	     'spline o1=0 d1=%g n1=%d' % (d,int(1.5+xmax/d)))
+	Flow(dipsfile,interfaces,'deriv scale=y')
+
+	Flow(modelfile,interfaces,
+	     '''
+	     unif2 d1=%g n1=%d v00=%s 
+	     ''' % (d,int(1.5+zmax/d),vstr))
+
+	Flow(raysfile,modelfile,
+		'''
+		rays2 yshot=3 zshot=1 nt=10000 dt=0.01 a0=50 amax=200 nr=1 escvar=n
+		''')
+
+#	Flow(huygens,modelfile,
+#		'''
+#		smooth repeat=3 rect1=5 rect2=3 |
+#		hwt2d xsou=0 zsou=3 nt=1000 d1=0.01 ng=1801 og=90 dg=1 rays=y
+#		''')
+
+	# plot the model
+	Plot(modelfile,'grey color=j scalebar=y label1=Depth unit1=km label2=Position unit2=km barlabel=Velocity barunit=km/s barreverse=y title=Model allpos=y')
+
+	# plot the ray
+	Plot(raysfile,'graph transp=y yreverse=y min1=0 max1=2 min2=0 max2=6 wantaxis=n wanttitle=n scalebar=y plotcol=7 plotfat=3')
+
+	# plot the ray
+#	Plot(huygens,'graph transp=n yreverse=y min1=0 max1=2 min2=0 max2=6 wantaxis=n wanttitle=n scalebar=y plotcol=7 plotfat=3')
+
+
+	# overlay model and ray
+	Result(overlay,[modelfile, raysfile],'Overlay')
+
+# Show error message if fail
+message = '''
+SCosncript Failed build
+'''
+
+atexit.register(print_build_failures,message)
+End()

experiments/modelUpdateStudy/errorReport.py

@@ -0,0 +1,54 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# errorReport.py  (Madagascar Recipe)
+#
+# Purpose: SCons error report function.
+# 
+# Important!: It should be called from a SConstruct 
+#
+# Site: http://www.dirackslounge.online
+# 
+# Version 1.0
+#
+# Programer: Rodolfo A. C. Neves (Dirack) 01/06/2020
+#
+# Email: rodolfo_profissional@hotmail.com
+#
+# License: GPL-3.0 <https://www.gnu.org/licenses/gpl-3.0.txt>.
+
+# Selfdoc string
+'''
+Python function to report errors in SCons
+
+If build fail report specific message. Import this module and add the following lineis to your SConstruct:
+
+import atexit
+from madagascarRecipes.errorReport import print_build_failures
+
+Call error message by the following:
+
+atexit.register(print_build_failures,message)
+
+It will call function print_build_failures to report message variable that should be a string with the error message to be reported
+'''
+
+if __name__ == "__main__":
+	print(__doc__)
+	exit()
+
+__author__="Rodolfo Dirack <rodolfo_profissional@hotmail.com>"
+__version__="1.0"
+
+def print_build_failures(message):
+	'''
+	Print error message when build fail
+	:param message: string, error message
+	'''
+	from SCons.Script import GetBuildFailures
+	for bf in GetBuildFailures():
+		print("BUILD ERROR REPORT".center(40,'*'))
+		print("FAILED BUILD: %s" % (bf.node))
+		print("MESSAGE REPORT".center(40,'*'))
+		print("%s" % (message))
+