Source code for pysit.vis.seismogram
import itertools
import time
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
__all__ = ['plot_seismogram']
[docs]def plot_seismogram(shot, axis=None, **kwargs):
if axis is None:
ax = plt.gca()
else:
ax = axis
domain = shot.sources.domain
data = shot.receivers.data
#np.array([r.data for r in shot.receiver_list])
plt.imshow(data, aspect='auto', **kwargs)
seismogram_tickers(shot)
plt.xlabel('X-Position')
plt.ylabel('Time (s)')
plt.show()
class ReceiverFormatterHelper(object):
def __init__(self, shot):
self.shot = shot
def __call__(self, grid_point, pos):
if ((grid_point >= 0) and (grid_point < len(self.shot.receivers.receiver_list))):
return '{0:.3}'.format(self.shot.receivers.receiver_list[int(np.floor(grid_point))].position[0])
else:
return ''
class TimeFormatterHelper(object):
def __init__(self, shot):
self.shot=shot
def __call__(self, grid_point, pos):
return '{0:.3}'.format(grid_point * self.shot.dt)
def seismogram_tickers(shot):
ax = plt.gca()
xformatter = mpl.ticker.FuncFormatter(ReceiverFormatterHelper(shot))
ax.xaxis.set_major_formatter(xformatter)
zformatter = mpl.ticker.FuncFormatter(TimeFormatterHelper(shot))
ax.yaxis.set_major_formatter(zformatter)
if __name__ == '__main__':
from pysit import *
from pysit.gallery import horizontal_reflector
# Setup
# Define Domain
pmlx = PML(0.1, 100)
pmlz = PML(0.1, 100)
x_config = (0.1, 1.0, pmlx, pmlx)
z_config = (0.1, 0.8, pmlz, pmlz)
d = RectangularDomain(x_config, z_config)
m = CartesianMesh(d, 90, 70)
# Generate true wave speed
C, C0, m, d = horizontal_reflector(m)
# Set up shots
zmin = d.z.lbound
zmax = d.z.rbound
zpos = zmin + (1./9.)*zmax
shots = equispaced_acquisition(m,
RickerWavelet(10.0),
sources=1,
source_depth=zpos,
source_kwargs={},
receivers='max',
receiver_depth=zpos,
receiver_kwargs={},
)
# Define and configure the wave solver
trange = (0.0,3.0)
solver = ConstantDensityAcousticWave(m,
# formulation='ode',
formulation='scalar',
model_parameters={'C': C},
spatial_accuracy_order=2,
trange=trange,
use_cpp_acceleration=True,
time_accuracy_order=4)
# Generate synthetic Seismic data
tt = time.time()
wavefields = []
base_model = solver.ModelParameters(m,{'C': C})
generate_seismic_data(shots, solver, base_model, wavefields=wavefields)
print('Data generation: {0}s'.format(time.time()-tt))
plot_seismogram(shots[0])