add examples and inversion funs

add examples and inversion funs

examples/seis_epic2d.py

@@ -0,0 +1,61 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Apr 29 09:19:19 2020
+
+@author: chens
+
+Seismic: 2D epicenter estimation assuming a homogeneous and flat Earth
+"""
+
+import numpy
+from geoist.pfm import giutils
+from geoist.inversion.geometry import Square
+import matplotlib.pyplot as plt
+from geoist.vis import giplt
+from geoist.inversion import ttime2d
+from geoist.inversion.epic2d import Homogeneous
+
+# Make a velocity model to calculate traveltimes
+area = (0, 10, 0, 10)
+vp, vs = 2, 1
+model = [Square(area, props={'vp': vp, 'vs': vs})]
+
+src = (5, 5)
+srcs = [src, src, src, src]
+recs = [(1, 2),(3,6),(4,7),(2,8)]
+
+#giutils.connect_points(src, rec_points)
+ptime = ttime2d.straight(model, 'vp', srcs, recs)
+stime = ttime2d.straight(model, 'vs', srcs, recs)
+# Calculate the residual time (S - P) with added noise
+traveltime, error = giutils.contaminate(stime - ptime, 0.05, percent=True,
+                                      return_stddev=True)
+solver = Homogeneous(traveltime, recs, vp, vs)
+
+initial = (1, 1)    
+estimate = solver.config('levmarq', initial=initial).fit().estimate_
+
+plt.figure(figsize=(10, 4))
+plt.subplot(1, 2, 1)
+plt.title('Epicenter + %d recording stations' % (len(recs)))
+plt.axis('scaled')
+giplt.points(src, '*y', label="True")
+giplt.points(recs, '^r', label="Stations")
+giplt.points(initial, '*b', label="Initial")
+giplt.points([estimate], '*g', label="Estimate")
+giplt.set_area(area)
+plt.legend(loc='lower right', shadow=True, numpoints=1, prop={'size': 12})
+plt.xlabel("X")
+plt.ylabel("Y")
+ax = plt.subplot(1, 2, 2)
+plt.title('Travel-time residuals + error bars')
+s = numpy.arange(len(traveltime)) + 1
+width = 0.3
+plt.bar(s - width, traveltime, width, color='g', label="Observed",
+        yerr=error)
+plt.bar(s, solver.predicted(), width, color='r', label="Predicted")
+ax.set_xticks(s)
+plt.legend(loc='upper right', shadow=True, prop={'size': 12})
+plt.xlabel("Station number")
+plt.ylabel("Travel-time residual")
+plt.show()

examples/seis_wavefd.py

@@ -0,0 +1,81 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Apr 29 09:45:57 2020
+
+@author: chens
+
+Seismic: 2D finite difference simulation of elastic P and SV wave propagation
+"""
+import numpy as np
+from matplotlib import animation  #from celluloid import Camera  #gif output
+import matplotlib.pyplot as plt
+from geoist.inversion import wavefd
+from geoist import DATA_PATH
+from pathlib import Path
+
+# Set the parameters of the finite difference grid
+shape = (200, 200)
+area = [0, 60000, 0, 60000]
+# Make a density and S wave velocity model
+density = 2400 * np.ones(shape)
+pvel = 6600
+svel = 3700
+mu = wavefd.lame_mu(svel, density)
+lamb = wavefd.lame_lamb(pvel, svel, density)
+
+# Make a wave source from a mexican hat wavelet that vibrates in the x and z
+# directions equaly
+sources = [[wavefd.MexHatSource(30000, 40000, area, shape, 10000, 1, delay=1)],
+           [wavefd.MexHatSource(30000, 40000, area, shape, 10000, 1, delay=1)]]
+
+# Get the iterator for the simulation
+dt = wavefd.maxdt(area, shape, pvel)
+duration = 20
+maxit = int(duration / dt)
+stations = [[55000, 0]]  # x, z coordinate of the seismometer
+snapshot = int(0.5 / dt)  # Plot a snapshot of the simulation every 0.5 seconds
+simulation = wavefd.elastic_psv(lamb, mu, density, area, dt, maxit, sources,
+                                stations, snapshot, padding=50, taper=0.01,
+                                xz2ps=True)
+
+# This part makes an animation using matplotlibs animation API
+fig = plt.figure(figsize=(12, 5))
+plt.subplot(2, 2, 2)
+plt.title('x component')
+xseismogram, = plt.plot([], [], '-k')
+plt.xlim(0, duration)
+plt.ylim(-10 ** (-3), 10 ** (-3))
+plt.subplot(2, 2, 4)
+plt.title('z component')
+zseismogram, = plt.plot([], [], '-k')
+plt.xlim(0, duration)
+plt.ylim(-10 ** (-3), 10 ** (-3))
+plt.subplot(1, 2, 1)
+# Start with everything zero and grab the plot so that it can be updated later
+wavefield = plt.imshow(np.zeros(shape), extent=area, vmin=-10 ** -6,
+                       vmax=10 ** -6, cmap=plt.cm.gray_r)
+plt.plot(stations[0][0], stations[0][1], '^k')
+plt.ylim(area[2:][::-1])
+plt.xlabel('x (km)')
+plt.ylabel('z (km)')
+times = np.linspace(0, maxit * dt, maxit)
+# This function updates the plot every few timesteps
+
+
+def animate(i):
+    """
+    Simulation will yield panels corresponding to P and S waves because
+    xz2ps=True
+    """
+    t, p, s, xcomp, zcomp = next(simulation)
+    plt.title('time: %0.1f s' % (times[t]))
+    wavefield.set_array((p + s)[::-1])
+    xseismogram.set_data(times[:t + 1], xcomp[0][:t + 1])
+    zseismogram.set_data(times[:t + 1], zcomp[0][:t + 1])
+    return wavefield, xseismogram, zseismogram
+
+
+anim = animation.FuncAnimation(fig, animate, frames= int(maxit / snapshot), interval=1)
+#animation.writers.list()
+anim.save(Path(DATA_PATH,'psv_wave.gif'),  writer = 'pillow', fps=20, dpi=200, bitrate=4000)
+plt.show()