import time
import numpy as np
import shutil #folder deletion
from scipy.interpolate import interp1d
import os
__all__ = ['generate_seismic_data', 'generate_seismic_data_from_file', 'generate_shot_data_time', 'generate_shot_data_frequency']
__docformat__ = "restructuredtext en"
[docs]def generate_seismic_data(shots, solver, model, verbose=False, frequencies=None, save_method=None, **kwargs):
"""Given a list of shots and a solver, generates seismic data.
Parameters
----------
shots : list of pysit.Shot
Collection of shots to be processed
solver : pysit.WaveSolver
Instance of wave solver to be used.
**kwargs : dict, optional
Optional arguments.
Notes
-----
`kwargs` may be used to specify `C0` and `wavefields` to `generate_shot_data`.
"""
if verbose:
print('Generating data...')
tt = time.time()
if solver.supports['equation_dynamics'] == "time":
for shot in shots:
generate_shot_data_time(shot, solver, model, verbose=verbose, **kwargs)
elif solver.supports['equation_dynamics'] == "frequency":
if frequencies is None:
raise TypeError('A frequency solver is passed, but no frequencies are given')
elif 'petsc' in kwargs and kwargs['petsc'] is not None:
# solve the Helmholtz operator for several rhs
generate_shot_data_frequency_list(shots, solver, model, frequencies, verbose=verbose, **kwargs)
else:
for shot in shots:
generate_shot_data_frequency(shot, solver, model, frequencies, verbose=verbose, **kwargs)
else:
raise TypeError("A time or frequency solver must be specified.")
if verbose:
data_tt = time.time() - tt
print('Data generation: {0}s'.format(data_tt))
print('Data generation: {0}s/shot'.format(data_tt/len(shots)))
if verbose:
print('Saving data...')
tt = time.time()
if solver.supports['equation_dynamics'] == "time":
# store shots in memory for a next simulation
if save_method is not None:
if save_method=='pickle':
try:
import pickle
except ImportError:
raise ImportError('cPickle is not installed please install it and try again')
newpath = r'./shots_time'
if os.path.exists(newpath):
shutil.rmtree(newpath)
if not os.path.exists(newpath):
os.makedirs(newpath)
for i in range(len(shots)):
fshot =newpath+'/shot_'+str(i+1)+'.save'
f = file(fshot, 'wb')
pickle.dump(shots[i].receivers.data,f, protocol=pickle.HIGHEST_PROTOCOL)
f.close()
fts =newpath+'/ts_'+str(i+1)+'.save'
f = file(fts, 'wb')
pickle.dump(shots[i].receivers.ts,f, protocol=pickle.HIGHEST_PROTOCOL)
f.close()
elif save_method=='savemat':
try:
import scipy.io as io
except ImportError:
raise ImportError('scipy.io is not installed please install it and try again')
newpath = r'./shots_time'
if os.path.exists(newpath):
shutil.rmtree(newpath)
if not os.path.exists(newpath):
os.makedirs(newpath)
for i in range(len(shots)):
fshot =newpath+'/shot_'+str(i+1)+'.mat'
io.savemat(fshot,mdict={fshot:shots[i].receivers.data})
fts =newpath+'/ts_'+str(i+1)+'.mat'
io.savemat(fts,mdict={fts:shots[i].receivers.ts})
elif save_method=='h5py':
try:
import h5py
except ImportError:
raise ImportError('h5py is not installed please install it and try again')
newpath = r'./shots_time'
if os.path.exists(newpath):
shutil.rmtree(newpath)
if not os.path.exists(newpath):
os.makedirs(newpath)
for i in range(len(shots)):
fshot =newpath+'/shot_'+str(i+1)+'.hdf5'
f = h5py.File(fshot,"w")
dts = f.create_dataset(fshot,shots[i].receivers.data.shape,shots[i].receivers.data.dtype)
dts[...] = shots[i].receivers.data
f.close()
fts =newpath+'/ts_'+str(i+1)+'.hdf5'
f = h5py.File(fts,'w')
dts = f.create_dataset(fts,shots[i].receivers.ts.shape,shots[i].receivers.ts.dtype)
dts[...] = shots[i].receivers.ts
f.close()
else:
raise TypeError('Unknown save_method')
if solver.supports['equation_dynamics'] == "frequency":
# store shots in memory for a next simulation
if save_method is not None:
if save_method=='pickle':
try:
import pickle
except ImportError:
raise ImportError('cPickle is not installed please install it and try again')
newpath = r'./shots_frequency'
if os.path.exists(newpath):
shutil.rmtree(newpath)
if not os.path.exists(newpath):
os.makedirs(newpath)
for i in range(len(shots)):
fshot =newpath+'/shot_'+str(i+1)+'.save'
f = file(fshot, 'wb')
pickle.dump(shots[i].receivers.data_dft,f, protocol=pickle.HIGHEST_PROTOCOL)
f.close()
elif save_method=='savemat':
try:
import scipy.io as io
except ImportError:
raise ImportError('scipy.io is not installed please install it and try again')
newpath = r'./shots_frequency'
if os.path.exists(newpath):
shutil.rmtree(newpath)
if not os.path.exists(newpath):
os.makedirs(newpath)
for i in range(len(shots)):
data_frequency = shots[i].receivers.data_dft
for nu in data_frequency:
fshot =newpath+'/shot_'+str(i+1)+'_nu_'+str(nu)+'.mat'
io.savemat(fshot,mdict={fshot:shots[i].receivers.data_dft[nu]})
elif save_method=='h5py':
raise NotImplementedError('h5py is not efficient with frequency data please use another data container')
else:
raise TypeError('Unknown save_method')
if verbose:
save_tt = time.time() - tt
print('Data saving: {0}s'.format(save_tt))
print('Data saving: {0}s/shot'.format(save_tt/len(shots)))
[docs]def generate_seismic_data_from_file(shots, solver, verbose=False, save_method=None, **kwargs):
if verbose:
print('Loading data...')
tt = time.time()
if save_method is not None:
if solver.supports['equation_dynamics'] == "time":
if save_method=='pickle':
try:
import pickle
except ImportError:
raise ImportError('cPickle is not installed please install it and try again')
path = r'./shots_time'
if os.path.exists(path):
try:
for i in range(len(shots)):
filename =path+'/ts_'+str(i+1)+'.save'
f = file(filename, 'rb')
ts = pickle.load(f)
f.close()
shots[i].receivers.clear_data(len(ts))
shots[i].receivers.ts = ts
filename =path+'/shot_'+str(i+1)+'.save'
f = file(filename, 'rb')
shots[i].receivers.data = pickle.load(f)
f.close()
shots[i].receivers.interpolator = interp1d(ts, np.zeros_like(shots[i].receivers.data),
axis=0, kind='linear', copy=False, bounds_error=False,
fill_value=0.0)
except IOError:
raise IOError('Files are corrupted or generated with different save_method argument')
else:
raise ImportError('There is no shot directory, please relaunch generate_seismic_data with save_method argument')
elif save_method=='savemat':
try:
import scipy.io as io
except ImportError:
raise ImportError('scipy.io is not installed please install it and try again')
path = r'./shots_time'
if os.path.exists(path):
try:
for i in range(len(shots)):
filename = path+'/ts_'+str(i+1)+'.mat'
ts = io.loadmat(filename)
ts = ts[filename]
ts = ts.flatten()
shots[i].receivers.clear_data(len(ts))
shots[i].receivers.ts = ts
filename = path+'/shot_'+str(i+1)+'.mat'
data = io.loadmat(filename)
data = data[filename]
shots[i].receivers.data = data
shots[i].receivers.interpolator = interp1d(ts, np.zeros_like(shots[i].receivers.data),
axis=0, kind='linear', copy=False, bounds_error=False,
fill_value=0.0)
except IOError:
raise IOError('Files are corrupted or generated with different save_method argument')
else:
raise ImportError('There is no shot directory, please relaunch generate_seismic_data with save_method argument')
elif save_method=='h5py':
try:
import h5py
except ImportError:
raise ImportError('h5py is not installed please install it and try again')
path = r'./shots_time'
if os.path.exists(path):
try:
for i in range(len(shots)):
filename =path+'/ts_'+str(i+1)+'.hdf5'
f = h5py.File(filename, 'r')
dts = f[filename]
ts = dts[...]
f.close()
shots[i].receivers.clear_data(len(ts))
shots[i].receivers.ts = ts
filename = path+'/shot_'+str(i+1)+'.hdf5'
f = h5py.File(filename, 'r')
dts = f[filename]
shots[i].receivers.data = dts[...]
f.close()
shots[i].receivers.interpolator = interp1d(ts, np.zeros_like(shots[i].receivers.data),
axis=0, kind='linear', copy=False, bounds_error=False,
fill_value=0.0)
except IOError:
raise IOError('Files are corrupted or generated with different save_method argument')
else:
raise ImportError('There is no shot directory, please relaunch generate_seismic_data with save_method argument')
else:
raise TypeError('Unknown save_method')
elif solver.supports['equation_dynamics'] == "frequency":
if save_method=='pickle':
try:
import pickle
except ImportError:
raise ImportError('cPickle is not installed please install it and try again')
path = r'./shots_frequency'
if os.path.exists(path):
try:
for i in range(len(shots)):
filename =path+'/shot_'+str(i+1)+'.save'
f = file(filename, 'rb')
shots[i].receivers.data_dft = pickle.load(f)
f.close()
except IOError:
raise IOError('Files are corrupted or generated with different save_method argument')
else:
raise ImportError('There is no shot directory, please relaunch generate_seismic_data with save_method argument')
elif save_method=='savemat':
try:
import scipy.io as io
except ImportError:
raise ImportError('scipy.io is not installed please install it and try again')
path = r'./shots_frequency'
if os.path.exists(path):
try:
if 'frequencies' in kwargs :
data_frequency = kwargs['frequencies']
for i in range(len(shots)):
shots[i].receivers.data_dft = dict()
for nu in data_frequency:
fshot =path+'/shot_'+str(i+1)+'_nu_'+str(nu)+'.mat'
data = io.loadmat(fshot)
data = data[fshot]
shots[i].receivers.data_dft[nu] = data
else :
raise AttributeError("""No frequencies specified in generate_seismic_data_from_file arguments""")
except IOError:
raise IOError('Files are corrupted or generated with different save_method argument')
else:
raise ImportError('There is no shot directory, please relaunch generate_seismic_data with save_method argument')
elif save_method=='h5py':
raise NotImplementedError('h5py is not efficient with frequency data please use another data container')
else:
raise AttributeError('No solver specified cannot load data')
else:
raise AttributeError('No save_method specified cannot load data')
if verbose:
load_tt = time.time() - tt
print('Data Loading: {0}s'.format(load_tt))
print('Data Loading: {0}s/shot'.format(load_tt/len(shots)))
[docs]def generate_shot_data_time(shot, solver, model, wavefields=None, wavefields_padded=None, verbose=False, **kwargs):
"""Given a shots and a solver, generates seismic data at the specified
receivers.
Parameters
----------
shots : list of pysit.Shot
Collection of shots to be processed
solver : pysit.WaveSolver
Instance of wave solver to be used.
model_parameters : solver.ModelParameters, optional
Wave equation parameters used for generating data.
wavefields : list, optional
List of wave states.
verbose : boolean
Verbosity flag.
Notes
-----
An empty list passed as `wavefields` will be populated with the state of the wave
field at each time index.
"""
solver.model_parameters = model
# Ensure that the receiver data is empty. And that interpolator is setup.
# shot.clear_data(solver.nsteps)
ts = solver.ts()
shot.reset_time_series(ts)
# Populate some stuff that will probably not exist soon anyway.
shot.dt = solver.dt
shot.trange = solver.trange
if solver.supports['equation_dynamics'] != "time":
raise TypeError('Solver must be a time solver to generate data.')
if(wavefields is not None):
wavefields[:] = []
if(wavefields_padded is not None):
wavefields_padded[:] = []
#Frequently used local references are faster than remote references
mesh = solver.mesh
dt = solver.dt
source = shot.sources
# Step k = 0
# p_0 is a zero array because if we assume the input signal is causal and we
# assume that the initial system (i.e., p_(-2) and p_(-1)) is uniformly
# zero, then the leapfrog scheme would compute that p_0 = 0 as well.
# This is stored in this solver specific data structure.
solver_data = solver.SolverData()
rhs_k = np.zeros(mesh.shape(include_bc=True))
rhs_kp1 = np.zeros(mesh.shape(include_bc=True))
# k is the t index. t = k*dt.
for k in range(solver.nsteps):
# print " Computing step {0}...".format(k)
uk = solver_data.k.primary_wavefield
# Extract the primary, non ghost or boundary nodes from the wavefield
# uk_bulk is a view so no new storage
# also, this is a bad way to do things. Somehow the BC padding needs to be better hidden. Maybe it needs to always be a part of the mesh?
uk_bulk = mesh.unpad_array(uk)
# Record the data at t_k
shot.receivers.sample_data_from_array(uk_bulk, k, **kwargs)
if(wavefields is not None):
wavefields.append(uk_bulk.copy())
if(wavefields_padded is not None):
wavefields_padded.append(uk.copy())
# When k is the nth step, the next time step is not needed, so save
# computation and break out early.
if(k == (solver.nsteps-1)): break
if k == 0:
rhs_k = mesh.pad_array(source.f(k*dt), out_array=rhs_k)
rhs_kp1 = mesh.pad_array(source.f((k+1)*dt), out_array=rhs_kp1)
else:
# shift time forward
rhs_k, rhs_kp1 = rhs_kp1, rhs_k
rhs_kp1 = mesh.pad_array(source.f((k+1)*dt), out_array=rhs_kp1)
# Given the state at k and k-1, compute the state at k+1
solver.time_step(solver_data, rhs_k, rhs_kp1)
# Don't know what data is needed for the solver, so the solver data
# handles advancing everything forward by one time step.
# k-1 <-- k, k <-- k+1, etc
solver_data.advance()
[docs]def generate_shot_data_frequency(shot, solver, model, frequencies, verbose=False, **kwargs):
"""most of this is copied from frequency_modeling.forward_model
Parameters
----------
shot : pysit.Shot
Gives the source signal approximation for the right hand side.
frequencies : list of 2-tuples
2-tuple, first element is the frequency to use, second element the weight.
return_parameters : list of {'wavefield', 'simdata', 'simdata_time', 'dWaveOp'}
Returns
-------
retval : dict
Dictionary whose keys are return_parameters that contains the specified data.
Notes
-----
* u is used as the target field universally. It could be velocity potential, it could be displacement, it could be pressure.
* uhat is used to generically refer to the DFT of u that is needed to compute the imaging condition.
"""
# Local references
solver.model_parameters = model
mesh = solver.mesh
source = shot.sources
# Sanitize the input
if not np.iterable(frequencies):
frequencies = [frequencies]
solver_data = solver.SolverData()
rhs = solver.WavefieldVector(mesh,dtype=solver.dtype)
for nu in frequencies:
rhs = solver.build_rhs(mesh.pad_array(source.f(nu=nu)), rhs_wavefieldvector=rhs)
solver.solve(solver_data, rhs, nu)
uhat = solver_data.k.primary_wavefield
# Record the data at frequency nu
shot.receivers.sample_data_from_array(mesh.unpad_array(uhat), nu=nu)
def generate_shot_data_frequency_list(shots, solver, model, frequencies, verbose=False, **kwargs):
"""most of this is copied from frequency_modeling.forward_model
Parameters
----------
shots : list of pysit.Shot
Gives the source signal approximation for the right hand side.
frequencies : list of 2-tuples
2-tuple, first element is the frequency to use, second element the weight.
return_parameters : list of {'wavefield', 'simdata', 'simdata_time', 'dWaveOp'}
Returns
-------
retval : dict
Dictionary whose keys are return_parameters that contains the specified data.
Notes
-----
* u is used as the target field universally. It could be velocity potential, it could be displacement, it could be pressure.
* uhat is used to generically refer to the DFT of u that is needed to compute the imaging condition.
"""
#number of shot varaible
nshot = len(shots)
# Local references
solver.model_parameters = model
mesh = solver.mesh
# Sanitize the input
if not np.iterable(frequencies):
frequencies = [frequencies]
rhs = solver.WavefieldVector(mesh,dtype=solver.dtype)
# decalare a list for the Right Hand Side
RHS = list()
for nu in frequencies:
del RHS[:]
for k in range(nshot):
shot = shots[k]
source = shot.sources
rhs = solver.build_rhs(mesh.pad_array(source.f(nu=nu)), rhs_wavefieldvector=rhs)
RHS.append(rhs.data.copy())
Uhat = solver.solve_petsc_uhat(solver, RHS, nu, **kwargs)
for k in range(nshot):
shot = shots[k]
uhat = Uhat[:,k]
# Giving the good dimension to the array
uhat = np.expand_dims(uhat, axis=1)
# Record the data at frequency nu
shot.receivers.sample_data_from_array(mesh.unpad_array(uhat), nu=nu)