marcel/pylot
3
0
Fork 0
Code Issues 7 Pull Requests 1 Releases Wiki Activity

[merge] resolved conflicts in utils due to two-sided coding

Browse Source
This commit is contained in:
Sebastian Wehling-Benatelli 2016-05-25 14:05:25 +02:00
commit f8db6b1d9f
6 changed files with 422 additions and 231 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

195
pylot/core/active/activeSeismoPick.py
View File

@ -3,6 +3,11 @@ import sys
import numpy as np import numpy as np
from pylot.core.active import seismicshot from pylot.core.active import seismicshot
from pylot.core.active.surveyUtils import cleanUp from pylot.core.active.surveyUtils import cleanUp
import copy_reg
import types
from pylot.core.util.utils import worker, _pickle_method
copy_reg.pickle(types.MethodType, _pickle_method)
class Survey(object): class Survey(object):
def __init__(self, path, sourcefile, receiverfile, useDefaultParas=False): def __init__(self, path, sourcefile, receiverfile, useDefaultParas=False):
@ -21,13 +26,13 @@ class Survey(object):
self._sourcefile = sourcefile self._sourcefile = sourcefile
self._obsdir = path self._obsdir = path
self._generateSurvey() self._generateSurvey()
self._initiateFilenames() self._initiate_fnames()
if useDefaultParas == True: if useDefaultParas == True:
self.setParametersForAllShots() self.setParametersForAllShots()
self._removeAllEmptyTraces() self._removeAllEmptyTraces()
self._updateShots() self._updateShots()
def _initiateFilenames(self): def _initiate_fnames(self):
for shot in self.data.values(): for shot in self.data.values():
shot.setRecfile(self.getPath() + self.getReceiverfile()) shot.setRecfile(self.getPath() + self.getReceiverfile())
shot.setSourcefile(self.getPath() + self.getSourcefile()) shot.setSourcefile(self.getPath() + self.getSourcefile())
@ -74,7 +79,7 @@ class Survey(object):
but were set in the input files. but were set in the input files.
''' '''
logfile = 'updateShots.out' logfile = 'updateShots.out'
count = 0; count = 0
countTraces = 0 countTraces = 0
for shot in self.data.values(): for shot in self.data.values():
del_traceIDs = shot.updateTraceList() del_traceIDs = shot.updateTraceList()
@ -93,15 +98,8 @@ class Survey(object):
"on removed traces." % (logfile)) "on removed traces." % (logfile))
outfile.close() outfile.close()
def setArtificialPick(self, traceID, pick): def setParametersForAllShots(self, cutwindow=(0, 0.2), tmovwind=0.3,
''' tsignal=0.03, tgap=0.0007):
Sets an artificial pick for a traceID of all shots in the survey object.
(Commonly used to generate a travel time t = 0 at the source origin)
'''
for shot in self.data.values():
shot.setPick(traceID, pick)
def setParametersForAllShots(self, cutwindow=(0, 0.2), tmovwind=0.3, tsignal=0.03, tgap=0.0007):
if (cutwindow == (0, 0.2) and tmovwind == 0.3 and if (cutwindow == (0, 0.2) and tmovwind == 0.3 and
tsignal == 0.03 and tgap == 0.0007): tsignal == 0.03 and tgap == 0.0007):
print ("Warning: Standard values used for " print ("Warning: Standard values used for "
@ -136,8 +134,10 @@ class Survey(object):
if not shot in diffs.keys(): if not shot in diffs.keys():
diffs[shot] = {} diffs[shot] = {}
for traceID in shot.getTraceIDlist(): for traceID in shot.getTraceIDlist():
if shot.getPickFlag(traceID) == 1 and shot.getManualPickFlag(traceID) == 1: if shot.getPickFlag(traceID) == 1 and shot.getManualPickFlag(
diffs[shot][traceID] = shot.getPick(traceID) - shot.getManualPick(traceID) traceID) == 1:
diffs[shot][traceID] = shot.getPick(
traceID) - shot.getManualPick(traceID)
return diffs return diffs
def plotDiffs(self): def plotDiffs(self):
@ -146,14 +146,15 @@ class Survey(object):
difference between automatic and manual pick. difference between automatic and manual pick.
''' '''
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
diffs = []; diffs = []
dists = []; dists = []
mpicks = []; mpicks = []
picks = [] picks = []
diffsDic = self.getDiffsFromManual() diffsDic = self.getDiffsFromManual()
for shot in self.data.values(): for shot in self.data.values():
for traceID in shot.getTraceIDlist(): for traceID in shot.getTraceIDlist():
if shot.getPickFlag(traceID) == 1 and shot.getManualPickFlag(traceID) == 1: if shot.getPickFlag(traceID) == 1 and shot.getManualPickFlag(
traceID) == 1:
dists.append(shot.getDistance(traceID)) dists.append(shot.getDistance(traceID))
mpicks.append(shot.getManualPick(traceID)) mpicks.append(shot.getManualPick(traceID))
picks.append(shot.getPick(traceID)) picks.append(shot.getPick(traceID))
@ -165,13 +166,16 @@ class Survey(object):
fig = plt.figure() fig = plt.figure()
ax = fig.add_subplot(111) ax = fig.add_subplot(111)
sc_a = ax.scatter(dists, picks, c='0.5', s=10, edgecolors='none', label=labela, alpha=0.3) sc_a = ax.scatter(dists, picks, c='0.5', s=10, edgecolors='none',
sc = ax.scatter(dists, mpicks, c=diffs, s=5, edgecolors='none', label=labelm) label=labela, alpha=0.3)
sc = ax.scatter(dists, mpicks, c=diffs, s=5, edgecolors='none',
label=labelm)
cbar = plt.colorbar(sc, fraction=0.05) cbar = plt.colorbar(sc, fraction=0.05)
cbar.set_label(labelm) cbar.set_label(labelm)
ax.set_xlabel('Distance [m]') ax.set_xlabel('Distance [m]')
ax.set_ylabel('Time [s]') ax.set_ylabel('Time [s]')
ax.text(0.5, 0.95, 'Plot of all MANUAL picks', transform=ax.transAxes, horizontalalignment='center') ax.text(0.5, 0.95, 'Plot of all MANUAL picks', transform=ax.transAxes,
horizontalalignment='center')
plt.legend() plt.legend()
def plotHist(self, nbins=20, ax=None): def plotHist(self, nbins=20, ax=None):
@ -186,14 +190,24 @@ class Survey(object):
ax = fig.add_subplot(111) ax = fig.add_subplot(111)
for shot in self.data.values(): for shot in self.data.values():
for traceID in shot.getTraceIDlist(): for traceID in shot.getTraceIDlist():
if shot.getPickFlag(traceID) == 1 and shot.getManualPickFlag(traceID) == 1: if shot.getPickFlag(traceID) == 1 and shot.getManualPickFlag(
traceID) == 1:
diffs.append(self.getDiffsFromManual()[shot][traceID]) diffs.append(self.getDiffsFromManual()[shot][traceID])
hist = plt.hist(diffs, nbins, histtype='step', normed=True, stacked=True) hist = plt.hist(diffs, nbins, histtype='step', normed=True,
plt.title('Histogram of the differences between automatic and manual pick') stacked=True)
plt.title(
'Histogram of the differences between automatic and manual pick')
plt.xlabel('Difference in time (auto - manual) [s]') plt.xlabel('Difference in time (auto - manual) [s]')
return diffs return diffs
def pickAllShots(self, vmin=333, vmax=5500, folm=0.6, HosAic='hos', aicwindow=(10, 0)): def pickShot(self, shTr):
shotnumber, traceID = shTr
shot = self.getShotForShotnumber(shotnumber)
traceID, pick = shot.pickTrace(traceID)
return shotnumber, traceID, pick
def pickAllShots(self, vmin=333, vmax=5500, folm=0.6, HosAic='hos',
aicwindow=(10, 0)):
''' '''
Automatically pick all traces of all shots of the survey. Automatically pick all traces of all shots of the survey.
@ -211,50 +225,58 @@ class Survey(object):
''' '''
from datetime import datetime from datetime import datetime
starttime = datetime.now() starttime = datetime.now()
count = 0; count = 0
tpicksum = starttime - starttime tpicksum = starttime - starttime
shTr = []
for shot in self.data.values(): for shot in self.data.values():
tstartpick = datetime.now(); tstartpick = datetime.now()
shot.setVmin(vmin)
shot.setVmax(vmax)
count += 1 count += 1
#shot.pickParallel(folm)
shot.setPickParameters(folm = folm, method = HosAic, aicwindow = aicwindow)
for traceID in shot.getTraceIDlist(): for traceID in shot.getTraceIDlist():
distance = shot.getDistance(traceID) # receive distance shTr.append((shot.getShotnumber(), traceID))
pickwin_used = shot.getCut() picks = worker(self.pickShot, shTr, async = True)
cutwindow = shot.getCut()
# for higher distances use a linear vmin/vmax to cut out late/early regions with high noise for shotnumber, traceID, pick in picks.get():
if distance > 5.: self.getShotForShotnumber(shotnumber).setPick(traceID, pick)
pwleft = distance / vmax ################## TEST
pwright = distance / vmin
if pwright > cutwindow[1]:
pwright = cutwindow[1]
pickwin_used = (pwleft, pwright)
shot.setPickwindow(traceID, pickwin_used) # tpicksum += (datetime.now() - tstartpick);
shot.pickTraces(traceID, folm, HosAic, aicwindow) # picker # tpick = tpicksum / count
# tremain = (tpick * (len(self.getShotDict()) - count))
# tend = datetime.now() + tremain
# progress = float(count) / float(len(self.getShotDict())) * 100
# self._update_progress(shot.getShotname(), tend, progress)
self.setSNR()
self.setEarllate()
print('\npickAllShots: Finished\n')
ntraces = self.countAllTraces()
pickedtraces = self.countAllPickedTraces()
print('Picked %s / %s traces (%d %%)\n'
% (pickedtraces, ntraces,
float(pickedtraces) / float(ntraces) * 100.))
def setSNR(self):
for shot in self.data.values():
for traceID in shot.getTraceIDlist():
shot.setSNR(traceID) shot.setSNR(traceID)
# if shot.getSNR(traceID)[0] < snrthreshold: # if shot.getSNR(traceID)[0] < snrthreshold:
if shot.getSNR(traceID)[0] < shot.getSNRthreshold(traceID): if shot.getSNR(traceID)[0] < shot.getSNRthreshold(traceID):
shot.removePick(traceID) shot.removePick(traceID)
def setEarllate(self):
for shot in self.data.values():
for traceID in shot.getTraceIDlist():
# set epp and lpp if SNR > 1 (else earllatepicker cant set values) # set epp and lpp if SNR > 1 (else earllatepicker cant set values)
if shot.getSNR(traceID)[0] > 1: if shot.getSNR(traceID)[0] > 1:
shot.setEarllatepick(traceID) shot.setEarllatepick(traceID)
tpicksum += (datetime.now() - tstartpick);
tpick = tpicksum / count
tremain = (tpick * (len(self.getShotDict()) - count))
tend = datetime.now() + tremain
progress = float(count) / float(len(self.getShotDict())) * 100
self._update_progress(shot.getShotname(), tend, progress)
print('\npickAllShots: Finished\n')
ntraces = self.countAllTraces()
pickedtraces = self.countAllPickedTraces()
print('Picked %s / %s traces (%d %%)\n'
% (pickedtraces, ntraces, float(pickedtraces) / float(ntraces) * 100.))
def cleanBySPE(self, maxSPE): def cleanBySPE(self, maxSPE):
''' '''
Sets all picks as invalid if they exceed a certain value of the symmetric pick error. Sets all picks as invalid if they exceed a certain value of the symmetric pick error.
@ -376,8 +398,11 @@ class Survey(object):
pickedTraces += 1 pickedTraces += 1
info_dict[shot.getShotnumber()] = {'numtraces': numtraces, info_dict[shot.getShotnumber()] = {'numtraces': numtraces,
'picked traces': [pickedTraces, 'picked traces': [pickedTraces,
'%2.2f %%' % (float(pickedTraces) / '%2.2f %%' % (
float(numtraces) * 100)], float(
pickedTraces) /
float(
numtraces) * 100)],
'mean SNR': np.mean(snrlist), 'mean SNR': np.mean(snrlist),
'mean distance': np.mean(dist)} 'mean distance': np.mean(dist)}
@ -391,10 +416,12 @@ class Survey(object):
if shot.getShotnumber() == shotnumber: if shot.getShotnumber() == shotnumber:
return shot return shot
def exportFMTOMO(self, directory='FMTOMO_export', sourcefile='input_sf.in', ttFileExtension='.tt'): def exportFMTOMO(self, directory='FMTOMO_export', sourcefile='input_sf.in',
ttFileExtension='.tt'):
''' '''
Exports all picks into a directory as travel time files readable by FMTOMO obsdata. Exports all picks into a directory as travel time files readable by FMTOMO obsdata.
''' '''
def getAngle(distance): def getAngle(distance):
PI = np.pi PI = np.pi
R = 6371. R = 6371.
@ -403,18 +430,20 @@ class Survey(object):
count = 0 count = 0
fmtomo_factor = 1000 # transforming [m/s] -> [km/s] fmtomo_factor = 1000 # transforming [m/s] -> [km/s]
LatAll = []; LatAll = []
LonAll = []; LonAll = []
DepthAll = [] DepthAll = []
srcfile = open(directory + '/' + sourcefile, 'w') srcfile = open(directory + '/' + sourcefile, 'w')
srcfile.writelines('%10s\n' % len(self.data)) # number of sources srcfile.writelines('%10s\n' % len(self.data)) # number of sources
for shotnumber in self.getShotlist(): for shotnumber in self.getShotlist():
shot = self.getShotForShotnumber(shotnumber) shot = self.getShotForShotnumber(shotnumber)
ttfilename = str(shotnumber) + ttFileExtension # filename of travel time file for this shot ttfilename = str(
shotnumber) + ttFileExtension # filename of travel time file for this shot
(x, y, z) = shot.getSrcLoc() # getSrcLoc returns (x, y, z) (x, y, z) = shot.getSrcLoc() # getSrcLoc returns (x, y, z)
srcfile.writelines('%10s %10s %10s\n' % (getAngle(y), getAngle(x), (-1) * z)) # transform to lat, lon, depth srcfile.writelines('%10s %10s %10s\n' % (
LatAll.append(getAngle(y)); getAngle(y), getAngle(x), (-1) * z)) # transform to lat, lon, depth
LonAll.append(getAngle(x)); LatAll.append(getAngle(y))
LonAll.append(getAngle(x))
DepthAll.append((-1) * z) DepthAll.append((-1) * z)
srcfile.writelines('%10s\n' % 1) srcfile.writelines('%10s\n' % 1)
srcfile.writelines('%10s %10s %10s\n' % (1, 1, ttfilename)) srcfile.writelines('%10s %10s %10s\n' % (1, 1, ttfilename))
@ -427,9 +456,10 @@ class Survey(object):
pick = shot.getPick(traceID) * fmtomo_factor pick = shot.getPick(traceID) * fmtomo_factor
delta = shot.getSymmetricPickError(traceID) * fmtomo_factor delta = shot.getSymmetricPickError(traceID) * fmtomo_factor
(x, y, z) = shot.getRecLoc(traceID) (x, y, z) = shot.getRecLoc(traceID)
ttfile.writelines('%20s %20s %20s %10s %10s\n' % (getAngle(y), getAngle(x), (-1) * z, pick, delta)) ttfile.writelines('%20s %20s %20s %10s %10s\n' % (
LatAll.append(getAngle(y)); getAngle(y), getAngle(x), (-1) * z, pick, delta))
LonAll.append(getAngle(x)); LatAll.append(getAngle(y))
LonAll.append(getAngle(x))
DepthAll.append((-1) * z) DepthAll.append((-1) * z)
count += 1 count += 1
ttfile.close() ttfile.close()
@ -489,19 +519,24 @@ class Survey(object):
if index <= figPerSubplot: if index <= figPerSubplot:
ax = fig.add_subplot(rows, columns, index) ax = fig.add_subplot(rows, columns, index)
if mode == '3d': if mode == '3d':
self.getShot(shotnumber).matshow(ax=ax, colorbar=False, annotations=True, legend=False) self.getShot(shotnumber).matshow(ax=ax, colorbar=False,
annotations=True,
legend=False)
elif mode == '2d': elif mode == '2d':
self.getShot(shotnumber).plot2dttc(ax) self.getShot(shotnumber).plot2dttc(ax)
self.getShot(shotnumber).plotmanual2dttc(ax) self.getShot(shotnumber).plotmanual2dttc(ax)
index += 1 index += 1
if index > figPerSubplot: if index > figPerSubplot:
fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0) fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0,
hspace=0)
fig = plt.figure() fig = plt.figure()
index = 1 index = 1
fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0) fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0,
hspace=0)
def plotAllPicks(self, plotRemoved=False, colorByVal='log10SNR', ax=None, cbar=None, refreshPlot=False): def plotAllPicks(self, plotRemoved=False, colorByVal='log10SNR', ax=None,
cbar=None, refreshPlot=False):
''' '''
Plots all picks over the distance between source and receiver. Returns (ax, region). Plots all picks over the distance between source and receiver. Returns (ax, region).
Picks can be checked and removed by using region class (pylot.core.active.surveyPlotTools.regions) Picks can be checked and removed by using region class (pylot.core.active.surveyPlotTools.regions)
@ -545,7 +580,8 @@ class Survey(object):
for shot in self.data.values(): for shot in self.data.values():
for traceID in shot.getTraceIDlist(): for traceID in shot.getTraceIDlist():
if plotRemoved == False: if plotRemoved == False:
if shot.getPickFlag(traceID) is not 0 or plotRemoved == True: if shot.getPickFlag(
traceID) is not 0 or plotRemoved == True:
dist.append(shot.getDistance(traceID)) dist.append(shot.getDistance(traceID))
pick.append(shot.getPick(traceID)) pick.append(shot.getPick(traceID))
snrlog.append(math.log10(shot.getSNR(traceID)[0])) snrlog.append(math.log10(shot.getSNR(traceID)[0]))
@ -557,12 +593,15 @@ class Survey(object):
'spe': spe} 'spe': spe}
self.color = color self.color = color
if refreshPlot is False: if refreshPlot is False:
ax, cbar = self.createPlot(dist, pick, color[colorByVal], label='%s' % colorByVal) ax, cbar = self.createPlot(dist, pick, color[colorByVal],
label='%s' % colorByVal)
region = regions(ax, cbar, self) region = regions(ax, cbar, self)
ax.legend() ax.legend()
return (ax, region) return (ax, region)
if refreshPlot is True: if refreshPlot is True:
ax, cbar = self.createPlot(dist, pick, color[colorByVal], label='%s' % colorByVal, ax=ax, cbar=cbar) ax, cbar = self.createPlot(dist, pick, color[colorByVal],
label='%s' % colorByVal, ax=ax,
cbar=cbar)
ax.legend() ax.legend()
return ax return ax
@ -577,23 +616,29 @@ class Survey(object):
print('Generating new plot...') print('Generating new plot...')
fig = plt.figure() fig = plt.figure()
ax = fig.add_subplot(111) ax = fig.add_subplot(111)
sc = ax.scatter(dist, pick, cmap=cm, c=inkByVal, s=5, edgecolors='none', label=label) sc = ax.scatter(dist, pick, cmap=cm, c=inkByVal, s=5,
edgecolors='none', label=label)
cbar = plt.colorbar(sc, fraction=0.05) cbar = plt.colorbar(sc, fraction=0.05)
cbar.set_label(label) cbar.set_label(label)
ax.set_xlabel('Distance [m]') ax.set_xlabel('Distance [m]')
ax.set_ylabel('Time [s]') ax.set_ylabel('Time [s]')
ax.text(0.5, 0.95, 'Plot of all picks', transform=ax.transAxes, horizontalalignment='center') ax.text(0.5, 0.95, 'Plot of all picks', transform=ax.transAxes,
horizontalalignment='center')
else: else:
sc = ax.scatter(dist, pick, cmap=cm, c=inkByVal, s=5, edgecolors='none', label=label) sc = ax.scatter(dist, pick, cmap=cm, c=inkByVal, s=5,
edgecolors='none', label=label)
cbar = plt.colorbar(sc, cax=cbar.ax) cbar = plt.colorbar(sc, cax=cbar.ax)
cbar.set_label(label) cbar.set_label(label)
ax.set_xlabel('Distance [m]') ax.set_xlabel('Distance [m]')
ax.set_ylabel('Time [s]') ax.set_ylabel('Time [s]')
ax.text(0.5, 0.95, 'Plot of all picks', transform=ax.transAxes, horizontalalignment='center') ax.text(0.5, 0.95, 'Plot of all picks', transform=ax.transAxes,
horizontalalignment='center')
return (ax, cbar) return (ax, cbar)
def _update_progress(self, shotname, tend, progress): def _update_progress(self, shotname, tend, progress):
sys.stdout.write('Working on shot %s. ETC is %02d:%02d:%02d [%2.2f %%]\r' % (shotname, sys.stdout.write(
'Working on shot %s. ETC is %02d:%02d:%02d [%2.2f %%]\r' % (
shotname,
tend.hour, tend.hour,
tend.minute, tend.minute,
tend.second, tend.second,

200
pylot/core/active/fmtomoUtils.py
View File

@ -6,23 +6,19 @@ import datetime
import numpy as np import numpy as np
class Tomo3d(object): class Tomo3d(object):
def __init__(self, nproc, iterations, citer = 0, overwrite = False): def __init__(self, fmtomodir, simuldir = 'fmtomo_simulation', citer = 0, overwrite = False):
''' '''
Class build from FMTOMO script tomo3d. Can be used to run several instances of FMM code in parallel. Class build from FMTOMO script tomo3d. Can be used to run several instances of FMM code in parallel.
:param: nproc, number of parallel processes
:type: integer
:param: iterations, number of iterations
:type: integer
:param: citer, current iteration (default = 0: start new model) :param: citer, current iteration (default = 0: start new model)
:type: integer :type: integer
''' '''
self.simuldir = simuldir
self.setCWD() self.setCWD()
self.buildFmtomodir(fmtomodir)
self.buildObsdata()
self.defParas() self.defParas()
self.nproc = nproc self.copyRef()
self.iter = iterations # number of iterations
self.citer = citer # current iteration self.citer = citer # current iteration
self.sources = self.readSrcFile() self.sources = self.readSrcFile()
self.traces = self.readTraces() self.traces = self.readTraces()
@ -33,22 +29,61 @@ class Tomo3d(object):
self.defFMMParas() self.defFMMParas()
self.defInvParas() self.defInvParas()
def buildFmtomodir(self, directory):
tomo_files = ['fm3d',
'frechgen',
'frechgen.in',
'invert3d',
'invert3d.in',
'mode_set.in',
'obsdata',
'obsdata.in',
'residuals',
'residuals.in',
'tomo3d',
'tomo3d.in']
for name in tomo_files:
filename = os.path.join(directory, name)
linkname = os.path.join(self.cwd, name)
os.system('ln -s %s %s'%(filename, linkname))
def buildObsdata(self):
os.system('obsdata')
os.system('mv sources.in sourcesref.in')
def defFMMParas(self): def defFMMParas(self):
'''
Initiates parameters for the forward calculation.
'''
# Name of fast marching program
self.fmm = '{0}/fm3d'.format(self.cwd) self.fmm = '{0}/fm3d'.format(self.cwd)
# Name of program calculating Frechet derivatives
self.frechgen = '{0}/frechgen'.format(self.cwd) self.frechgen = '{0}/frechgen'.format(self.cwd)
# Name of current velocity/inversion grid
self.cvg = 'vgrids.in' self.cvg = 'vgrids.in'
# Name of current interfaces grid
self.cig = 'interfaces.in' self.cig = 'interfaces.in'
# Name of file containing current source locations
self.csl = 'sources.in' self.csl = 'sources.in'
# Name of file containing propagation grid
self.pg = 'propgrid.in' self.pg = 'propgrid.in'
# Name of file containing receiver coordinates
self.rec = 'receivers.in' self.rec = 'receivers.in'
self.frech = 'frechet.in' self.frech = 'frechet.in'
self.frechout = 'frechet.dat' self.frechout = 'frechet.dat'
# Name of file containing measured data
self.ot = 'otimes.dat' self.ot = 'otimes.dat'
# Name of file containing output velocity information
self.ttim = 'arrivals.dat' self.ttim = 'arrivals.dat'
self.mode = 'mode_set.in' self.mode = 'mode_set.in'
# Name of temporary folders created for each process
self.folder = '.proc_' self.folder = '.proc_'
def defInvParas(self): def defInvParas(self):
'''
Initiates inversion parameters for FMTOMO.
'''
# Name of program for performing inversion # Name of program for performing inversion
self.inv = '{0}/invert3d'.format(self.cwd) self.inv = '{0}/invert3d'.format(self.cwd)
# Name of file containing current model traveltimes # Name of file containing current model traveltimes
@ -67,19 +102,42 @@ class Tomo3d(object):
self.resout = 'residuals.dat' self.resout = 'residuals.dat'
def copyRef(self): def copyRef(self):
# Copies reference grids to used grids (e.g. sourcesref.in to sources.in) '''
Copies reference grids to used grids (e.g. sourcesref.in to sources.in)
'''
os.system('cp %s %s'%(self.ivg, self.cvg)) os.system('cp %s %s'%(self.ivg, self.cvg))
os.system('cp %s %s'%(self.iig, self.cig)) os.system('cp %s %s'%(self.iig, self.cig))
os.system('cp %s %s'%(self.isl, self.csl)) os.system('cp %s %s'%(self.isl, self.csl))
def setCWD(self): def setCWD(self, directory = None):
self.cwd = subprocess.check_output(['pwd'])[0:-1] '''
print('Working directory is pwd: %s'%self.cwd) Set working directory containing all necessary files.
Default: pwd
'''
if directory == None:
directory = os.path.join(os.getcwd(), self.simuldir)
os.chdir(directory)
self.cwd = directory
print('Working directory is: %s'%self.cwd)
def runFrech(self): def runFrech(self):
os.system(self.frechgen) os.system(self.frechgen)
def runTOMO3D(self): def runTOMO3D(self, nproc, iterations):
'''
Starts up the FMTOMO code for the set number of iterations on nproc parallel processes.
:param: nproc, number of parallel processes
:type: integer
:param: iterations, number of iterations
:type: integer
'''
self.nproc = nproc
self.iter = iterations # number of iterations
starttime = datetime.datetime.now() starttime = datetime.datetime.now()
print('Starting TOMO3D on %s parallel processes for %s iteration(s).' print('Starting TOMO3D on %s parallel processes for %s iteration(s).'
%(self.nproc, self.iter)) %(self.nproc, self.iter))
@ -104,6 +162,10 @@ class Tomo3d(object):
print('runTOMO3D: Finished %s iterations after %s.'%(self.iter, tdelta)) print('runTOMO3D: Finished %s iterations after %s.'%(self.iter, tdelta))
def runFmm(self, directory, logfile, processes): def runFmm(self, directory, logfile, processes):
'''
Calls an instance of the FMM code in the process directory.
Requires a list of all active processes and returns an updated list.
'''
os.chdir(directory) os.chdir(directory)
fout = open(logfile, 'w') fout = open(logfile, 'w')
processes.append(subprocess.Popen(self.fmm, stdout = fout)) processes.append(subprocess.Popen(self.fmm, stdout = fout))
@ -112,6 +174,9 @@ class Tomo3d(object):
return processes return processes
def startForward(self, logdir): def startForward(self, logdir):
'''
Runs an instance of the FMM code in the process directory.
'''
self._printLine() self._printLine()
print('Starting forward simulation for iteration %s.'%(self.citer)) print('Starting forward simulation for iteration %s.'%(self.citer))
@ -128,8 +193,8 @@ class Tomo3d(object):
logfn = 'fm3dlog_' + str(procID) + '.out' logfn = 'fm3dlog_' + str(procID) + '.out'
log_out = os.path.join(logdir, logfn) log_out = os.path.join(logdir, logfn)
self.writeSrcFile(procID, directory) self.writeSrcFile(procID)
self.writeTracesFile(procID, directory) self.writeTracesFile(procID)
os.system('cp {cvg} {cig} {mode} {pg} {frechin} {dest}' os.system('cp {cvg} {cig} {mode} {pg} {frechin} {dest}'
.format(cvg=self.cvg, cig=self.cig, frechin=self.frech, .format(cvg=self.cvg, cig=self.cig, frechin=self.frech,
mode=self.mode, pg=self.pg, dest=directory)) mode=self.mode, pg=self.pg, dest=directory))
@ -149,10 +214,16 @@ class Tomo3d(object):
print('Finished Forward calculation after %s'%tdelta) print('Finished Forward calculation after %s'%tdelta)
def startInversion(self): def startInversion(self):
'''
Simply calls the inversion program.
'''
print('Calling %s...'%self.inv) print('Calling %s...'%self.inv)
os.system(self.inv) os.system(self.inv)
def calcRes(self): def calcRes(self):
'''
Calls residual calculation program.
'''
resout = os.path.join(self.cwd, self.resout) resout = os.path.join(self.cwd, self.resout)
if self.citer == 0: if self.citer == 0:
os.system('%s > %s'%(self.resid, resout)) os.system('%s > %s'%(self.resid, resout))
@ -185,11 +256,17 @@ class Tomo3d(object):
raise RuntimeError('Could not create directory: %s'%directory) raise RuntimeError('Could not create directory: %s'%directory)
def makeDirectories(self): def makeDirectories(self):
'''
Makes temporary directories for all processes.
'''
for procID in range(1, self.nproc + 1): for procID in range(1, self.nproc + 1):
directory = self.getProcDir(procID) directory = self.getProcDir(procID)
self.makeDir(directory) self.makeDir(directory)
def makeInvIterDir(self): def makeInvIterDir(self):
'''
Makes directories for each iteration step for the output.
'''
invIterDir = self.cwd + '/it_%s'%(self.citer) invIterDir = self.cwd + '/it_%s'%(self.citer)
err = os.system('mkdir %s'%invIterDir) err = os.system('mkdir %s'%invIterDir)
if err is 256: if err is 256:
@ -200,12 +277,18 @@ class Tomo3d(object):
self.cInvIterDir = invIterDir self.cInvIterDir = invIterDir
def clearDir(self, directory): def clearDir(self, directory):
print('Wiping directory %s...'%directory) '''
Wipes a certain directory.
'''
#print('Wiping directory %s...'%directory)
for filename in os.listdir(directory): for filename in os.listdir(directory):
filenp = os.path.join(directory, filename) filenp = os.path.join(directory, filename)
os.remove(filenp) os.remove(filenp)
def clearDirectories(self): def clearDirectories(self):
'''
Wipes all generated temporary directories.
'''
for directory in self.directories: for directory in self.directories:
self.clearDir(directory) self.clearDir(directory)
@ -214,14 +297,24 @@ class Tomo3d(object):
return os.rmdir(directory) return os.rmdir(directory)
def removeDirectories(self): def removeDirectories(self):
'''
Removes all generated temporary directories.
'''
for directory in self.directories: for directory in self.directories:
self.rmDir(directory) self.rmDir(directory)
self.directories = [] self.directories = []
def getProcDir(self, procID): def getProcDir(self, procID):
'''
Returns the temporary directory for a certain process
with procID = process number.
'''
return os.path.join(self.cwd, self.folder) + str(procID) return os.path.join(self.cwd, self.folder) + str(procID)
def getTraceIDs4Sources(self, sourceIDs): def getTraceIDs4Sources(self, sourceIDs):
'''
Returns corresponding trace IDs for a set of given source IDs.
'''
traceIDs = [] traceIDs = []
for traceID in self.traces.keys(): for traceID in self.traces.keys():
if self.traces[traceID]['source'] in sourceIDs: if self.traces[traceID]['source'] in sourceIDs:
@ -229,6 +322,9 @@ class Tomo3d(object):
return traceIDs return traceIDs
def getTraceIDs4Source(self, sourceID): def getTraceIDs4Source(self, sourceID):
'''
Returns corresponding trace IDs for a source ID.
'''
traceIDs = [] traceIDs = []
for traceID in self.traces.keys(): for traceID in self.traces.keys():
if self.traces[traceID]['source'] == sourceID: if self.traces[traceID]['source'] == sourceID:
@ -236,22 +332,38 @@ class Tomo3d(object):
return traceIDs return traceIDs
def copyArrivals(self, target = None): def copyArrivals(self, target = None):
'''
Copies the FMM output file (self.ttim) to a specific target file.
Default target is self.mtrav (model travel times).
'''
if target == None: if target == None:
target = os.path.join(self.cwd, self.mtrav) target = os.path.join(self.cwd, self.mtrav)
os.system('cp %s %s'%(os.path.join( os.system('cp %s %s'%(os.path.join(
self.cInvIterDir, self.ttim), target)) self.cInvIterDir, self.ttim), target))
def saveVgrid(self): def saveVgrid(self):
vgpath = os.path.join(self.cwd, 'vgrids.in') '''
Saves the current velocity grid for the current iteration step.
'''
vgpath = os.path.join(self.cwd, self.cvg)
os.system('cp %s %s'%(vgpath, self.cInvIterDir)) os.system('cp %s %s'%(vgpath, self.cInvIterDir))
def calcSrcPerKernel(self): def calcSrcPerKernel(self):
'''
(Equally) distributes all sources depending on the number of processes (kernels).
Returns two integer values.
First: minimum number of sources for each process
Second: remaining sources (always less than number of processes)
'''
nsrc = self.readNsrc() nsrc = self.readNsrc()
if self.nproc > nsrc: if self.nproc > nsrc:
raise ValueError('Number of spawned processes higher than number of sources') raise ValueError('Number of spawned processes higher than number of sources')
return nsrc/self.nproc, nsrc%self.nproc return nsrc/self.nproc, nsrc%self.nproc
def srcIDs4Kernel(self, procID): def srcIDs4Kernel(self, procID):
'''
Calculates and returns all source IDs for a given process ID.
'''
proc = procID - 1 proc = procID - 1
nsrc = self.readNsrc() nsrc = self.readNsrc()
srcPK, remain = self.calcSrcPerKernel() srcPK, remain = self.calcSrcPerKernel()
@ -274,12 +386,18 @@ class Tomo3d(object):
return nsrc return nsrc
def readNtraces(self): def readNtraces(self):
'''
Reads the total number of traces from self.rec header.
'''
recfile = open(self.rec, 'r') recfile = open(self.rec, 'r')
nrec = int(recfile.readline()) nrec = int(recfile.readline())
recfile.close() recfile.close()
return nrec return nrec
def readSrcFile(self): def readSrcFile(self):
'''
Reads the whole sourcefile and returns structured information in a dictionary.
'''
nsrc = self.readNsrc() nsrc = self.readNsrc()
srcfile = open(self.csl, 'r') srcfile = open(self.csl, 'r')
@ -309,6 +427,10 @@ class Tomo3d(object):
return sources return sources
def readTraces(self): def readTraces(self):
'''
Reads the receiver input file and returns the information
in a structured dictionary.
'''
recfile = open(self.rec, 'r') recfile = open(self.rec, 'r')
ntraces = self.readNtraces() ntraces = self.readNtraces()
@ -330,8 +452,13 @@ class Tomo3d(object):
return traces return traces
def readArrivals(self, procID): def readArrivals(self, procID):
'''
Reads the arrival times from a temporary process directory,
changes local to global sourceIDs and traceIDs and returns
a list of arrival times.
'''
directory = self.getProcDir(procID) directory = self.getProcDir(procID)
arrfile = open(directory + '/arrivals.dat', 'r') arrfile = open(os.path.join(directory, self.ttim), 'r')
sourceIDs = self.srcIDs4Kernel(procID) sourceIDs = self.srcIDs4Kernel(procID)
arrivals = [] arrivals = []
@ -347,6 +474,10 @@ class Tomo3d(object):
return arrivals return arrivals
def readRays(self, procID): def readRays(self, procID):
'''
Reads rays output from a temporary process directory and returns
the information in a structured dictionary.
'''
directory = self.getProcDir(procID) directory = self.getProcDir(procID)
raysfile = open(directory + '/rays.dat', 'r') raysfile = open(directory + '/rays.dat', 'r')
sourceIDs = self.srcIDs4Kernel(procID) sourceIDs = self.srcIDs4Kernel(procID)
@ -385,8 +516,12 @@ class Tomo3d(object):
} }
return rays return rays
def writeSrcFile(self, procID, directory): def writeSrcFile(self, procID):
srcfile = open('%s/sources.in'%directory, 'w') '''
Writes a source input file for a process with ID = procID.
'''
directory = self.getProcDir(procID)
srcfile = open(os.path.join(directory, self.csl), 'w')
sourceIDs = self.srcIDs4Kernel(procID) sourceIDs = self.srcIDs4Kernel(procID)
srcfile.write('%s\n'%len(sourceIDs)) srcfile.write('%s\n'%len(sourceIDs))
@ -401,7 +536,11 @@ class Tomo3d(object):
srcfile.write('%s %s\n'%(int(interactions[0]), int(interactions[1]))) srcfile.write('%s %s\n'%(int(interactions[0]), int(interactions[1])))
srcfile.write('%s\n'%source['veltype']) srcfile.write('%s\n'%source['veltype'])
def writeTracesFile(self, procID, directory): def writeTracesFile(self, procID):
'''
Writes a receiver input file for a process with ID = procID.
'''
directory = self.getProcDir(procID)
recfile = open('%s/receivers.in'%directory, 'w') recfile = open('%s/receivers.in'%directory, 'w')
sourceIDs = self.srcIDs4Kernel(procID) sourceIDs = self.srcIDs4Kernel(procID)
traceIDs = self.getTraceIDs4Sources(sourceIDs) traceIDs = self.getTraceIDs4Sources(sourceIDs)
@ -417,9 +556,12 @@ class Tomo3d(object):
recfile.write('%s\n'%trace['path']) recfile.write('%s\n'%trace['path'])
def mergeArrivals(self, directory): def mergeArrivals(self, directory):
'''
Merges the arrival times for all processes to self.cInvIterDir.
'''
arrfn = os.path.join(directory, self.ttim) arrfn = os.path.join(directory, self.ttim)
arrivalsOut = open(arrfn, 'w') arrivalsOut = open(arrfn, 'w')
print('Merging arrivals.dat...') print('Merging %s...'%self.ttim)
for procID in range(1, self.nproc + 1): for procID in range(1, self.nproc + 1):
arrivals = self.readArrivals(procID) arrivals = self.readArrivals(procID)
for line in arrivals: for line in arrivals:
@ -428,6 +570,9 @@ class Tomo3d(object):
os.system('ln -fs %s %s'%(arrfn, os.path.join(self.cwd, self.ttim))) os.system('ln -fs %s %s'%(arrfn, os.path.join(self.cwd, self.ttim)))
def mergeRays(self, directory): def mergeRays(self, directory):
'''
Merges the ray paths for all processes to self.cInvIterDir.
'''
print('Merging rays.dat...') print('Merging rays.dat...')
with open(directory + '/rays.dat', 'w') as outfile: with open(directory + '/rays.dat', 'w') as outfile:
for procID in range(1, self.nproc + 1): for procID in range(1, self.nproc + 1):
@ -448,9 +593,11 @@ class Tomo3d(object):
outfile.writelines(raysec['raypoints']) outfile.writelines(raysec['raypoints'])
def mergeFrechet(self, directory): def mergeFrechet(self, directory):
print('Merging frechet.dat...') '''
Merges the frechet derivatives for all processes to self.cInvIterDir.
'''
frechfnout = os.path.join(directory, self.frechout) frechfnout = os.path.join(directory, self.frechout)
print('Merging %s...'%self.frechout)
with open(frechfnout, 'w') as outfile: with open(frechfnout, 'w') as outfile:
for procID in range(1, self.nproc + 1): for procID in range(1, self.nproc + 1):
filename = os.path.join(self.getProcDir(procID), self.frechout) filename = os.path.join(self.getProcDir(procID), self.frechout)
@ -465,6 +612,9 @@ class Tomo3d(object):
os.system('ln -fs %s %s'%(frechfnout, os.path.join(self.cwd, self.frechout))) os.system('ln -fs %s %s'%(frechfnout, os.path.join(self.cwd, self.frechout)))
def mergeOutput(self, directory): def mergeOutput(self, directory):
'''
Calls self.mergeArrivals, self.mergeFrechet and self.mergeRays.
'''
self.mergeArrivals(directory) self.mergeArrivals(directory)
self.mergeFrechet(directory) self.mergeFrechet(directory)
self.mergeRays(directory) self.mergeRays(directory)

8
pylot/core/active/seismicArrayPreparation.py
View File

@ -401,7 +401,7 @@ class SeisArray(object):
customgrid='mygrid.in', writeVTK=True): customgrid='mygrid.in', writeVTK=True):
''' '''
Generate FMTOMO input files from the SeisArray dimensions. Generate FMTOMO input files from the SeisArray dimensions.
Generates: vgrids.in, interfaces.in, propgrid.in Generates: vgridsref.in, interfacesref.in, propgrid.in
:param: nPointsPropgrid, number of points in each direction of the propagation grid (z, y, x) :param: nPointsPropgrid, number of points in each direction of the propagation grid (z, y, x)
:type: tuple :type: tuple
@ -473,10 +473,10 @@ class SeisArray(object):
outfile.close() outfile.close()
def generateInterfaces(self, nTheta, nPhi, depthmax, cushionfactor=0.1, def generateInterfaces(self, nTheta, nPhi, depthmax, cushionfactor=0.1,
outfilename='interfaces.in', method='linear', outfilename='interfacesref.in', method='linear',
returnInterfaces=False): returnInterfaces=False):
''' '''
Create an interfaces.in file for FMTOMO from the SeisArray boundaries. Create an interfacesref.in file for FMTOMO from the SeisArray boundaries.
:param: nTheta, number of points in Theta :param: nTheta, number of points in Theta
type: int type: int
@ -610,7 +610,7 @@ class SeisArray(object):
def generateVgrid(self, nTheta, nPhi, nR, Rbt, thetaSN=None, def generateVgrid(self, nTheta, nPhi, nR, Rbt, thetaSN=None,
phiWE=None, cushionfactor=0.1, phiWE=None, cushionfactor=0.1,
outfilename='vgrids.in', method='linear', outfilename='vgridsref.in', method='linear',
infilename='mygrid.in', returnTopo=False): infilename='mygrid.in', returnTopo=False):
''' '''
Generate a velocity grid for fmtomo regarding topography with a linear gradient starting at the topography with 0.34 [km/s]. Generate a velocity grid for fmtomo regarding topography with a linear gradient starting at the topography with 0.34 [km/s].

151
pylot/core/active/seismicshot.py
View File

@ -11,10 +11,15 @@ from pylot.core.pick.charfuns import AICcf
from pylot.core.pick.utils import getSNR from pylot.core.pick.utils import getSNR
from pylot.core.pick.utils import earllatepicker from pylot.core.pick.utils import earllatepicker
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import warnings
import copy_reg
import types
from pylot.core.util.utils import worker, _pickle_method
copy_reg.pickle(types.MethodType, _pickle_method)
plt.interactive('True') plt.interactive('True')
class SeismicShot(object): class SeismicShot(object):
''' '''
SuperClass for a seismic shot object. SuperClass for a seismic shot object.
@ -64,11 +69,6 @@ class SeismicShot(object):
if traceID == trace.stats.channel: if traceID == trace.stats.channel:
self.traces.remove(trace) self.traces.remove(trace)
# for traceID in TraceIDs:
# traces = [trace for trace in self.traces if int(trace.stats.channel) == traceID]
# if len(traces) is not 1:
# self.traces.remove(trace)
def updateTraceList(self): def updateTraceList(self):
''' '''
Looks for empty traces, returns a list of deleted traceIDs. Looks for empty traces, returns a list of deleted traceIDs.
@ -83,6 +83,12 @@ class SeismicShot(object):
def setParameters(self, name, value): def setParameters(self, name, value):
self.paras[name] = value self.paras[name] = value
def setVmin(self, vmin):
self.setParameters('vmin', vmin)
def setVmax(self, vmax):
self.setParameters('vmax', vmax)
def setCut(self, cut): def setCut(self, cut):
self.setParameters('cut', cut) self.setParameters('cut', cut)
@ -107,6 +113,43 @@ class SeismicShot(object):
def setSourcefile(self, sourcefile): def setSourcefile(self, sourcefile):
self.setParameters('sourcefile', sourcefile) self.setParameters('sourcefile', sourcefile)
def setMethod(self, method):
self.setParameters('method', method)
def setAicwindow(self, aicwindow):
self.setParameters('aicwindow', aicwindow)
def setFolm(self, folm):
self.setParameters('folm', folm)
def setDynPickwindow(self, traceID, cutdist = 5.):
distance = self.getDistance(traceID) # receive distance
vmin = self.getVmin()
vmax = self.getVmax()
pickwin_used = self.getCut()
cutwindow = self.getCut()
# for higher distances use a linear vmin/vmax to cut out late/early regions with high noise
if distance > cutdist:
pwleft = distance / vmax
pwright = distance / vmin
if pwright > cutwindow[1]:
pwright = cutwindow[1]
pickwin_used = (pwleft, pwright)
self.setPickwindow(traceID, pickwin_used)
def getMethod(self):
return self.paras['method']
def getAicwindow(self):
return self.paras['aicwindow']
def getFolm(self):
return self.paras['folm']
def getParas(self): def getParas(self):
return self.paras return self.paras
@ -137,6 +180,12 @@ class SeismicShot(object):
def getSourcefile(self): def getSourcefile(self):
return self.paras['sourcefile'] return self.paras['sourcefile']
def getVmin(self):
return self.paras['vmin']
def getVmax(self):
return self.paras['vmax']
def getManualPick(self, traceID): def getManualPick(self, traceID):
if not self.getManualPickFlag(traceID) == 0: if not self.getManualPickFlag(traceID) == 0:
return self.manualpicks[traceID]['mpp'] return self.manualpicks[traceID]['mpp']
@ -151,7 +200,6 @@ class SeismicShot(object):
if not self.getPickFlag(traceID) == 0: if not self.getPickFlag(traceID) == 0:
return self.picks[traceID]['mpp'] return self.picks[traceID]['mpp']
if returnRemoved == True: if returnRemoved == True:
# print('getPick: Returned removed pick for shot %d, traceID %d' %(self.getShotnumber(), traceID))
return self.picks[traceID]['mpp'] return self.picks[traceID]['mpp']
def getPickIncludeRemoved(self, traceID): def getPickIncludeRemoved(self, traceID):
@ -279,12 +327,28 @@ class SeismicShot(object):
if len(traces) == 1: if len(traces) == 1:
return Stream(traces) return Stream(traces)
self.setPick(traceID, None) self.setPick(traceID, None)
print 'Warning: ambigious or empty traceID: %s' % traceID warnings.warn('ambigious or empty traceID: %s' % traceID)
# raise ValueError('ambigious or empty traceID: %s' % traceID) def setPickParameters(self, folm, method = 'hos', aicwindow = (10, 0)):
self.setFolm(folm)
self.setMethod(method)
self.setAicwindow(aicwindow)
def pickTraces(self, traceID, windowsize, folm, HosAic='hos'): ########## input variables ########## # def pickParallel(self):
# LOCALMAX NOT IMPLEMENTED! # traceIDs = self.getTraceIDlist()
# picks = []
# #picks = worker(self.pickTrace, traceIDs)
# # for traceID, pick in picks:
# # self.setPick(traceID, pick)
# for traceID in traceIDs:
# trID, pick = self.pickTrace(traceID)
# picks.append([trID, pick])
# #self.setPick(traceID, pick)
# return picks
def pickTrace(self, traceID):
''' '''
Intitiate picking for a trace. Intitiate picking for a trace.
@ -309,17 +373,17 @@ class SeismicShot(object):
:param: HosAic, get hos or aic pick (can be 'hos'(default) or 'aic') :param: HosAic, get hos or aic pick (can be 'hos'(default) or 'aic')
:type: 'string' :type: 'string'
''' '''
self.setDynPickwindow(traceID)
hoscf = self.getHOScf(traceID) ### determination of both, HOS and AIC (need to change threshold-picker) ### hoscf = self.getHOScf(traceID) ### determination of both, HOS and AIC (need to change threshold-picker) ###
aiccf = self.getAICcf(traceID) aiccf = self.getAICcf(traceID)
self.folm = folm
self.timeArray[traceID] = hoscf.getTimeArray() self.timeArray[traceID] = hoscf.getTimeArray()
aiccftime, hoscftime = self.threshold(hoscf, aiccf, windowsize, self.getPickwindow(traceID), folm) aiccftime, hoscftime = self.threshold(hoscf, aiccf, self.getAicwindow(), self.getPickwindow(traceID), self.getFolm())
setHosAic = {'hos': hoscftime, setHosAic = {'hos': hoscftime,
'aic': aiccftime} 'aic': aiccftime}
self.setPick(traceID, setHosAic[HosAic]) return traceID, setHosAic[self.getMethod()]
def setEarllatepick(self, traceID, nfac=1.5): def setEarllatepick(self, traceID, nfac=1.5):
tgap = self.getTgap() tgap = self.getTgap()
@ -335,11 +399,6 @@ class SeismicShot(object):
if self.picks[traceID]['epp'] < 0: if self.picks[traceID]['epp'] < 0:
self.picks[traceID]['epp'] self.picks[traceID]['epp']
# print('setEarllatepick: Set epp to 0 because it was < 0')
# TEST OF 1/2 PICKERROR
# self.picks[traceID]['spe'] *= 0.5
# TEST OF 1/2 PICKERROR
def threshold(self, hoscf, aiccf, windowsize, pickwindow, folm): def threshold(self, hoscf, aiccf, windowsize, pickwindow, folm):
''' '''
@ -368,12 +427,8 @@ class SeismicShot(object):
leftb = int(pickwindow[0] / self.getCut()[1] * len(hoscflist)) leftb = int(pickwindow[0] / self.getCut()[1] * len(hoscflist))
rightb = int(pickwindow[1] / self.getCut()[1] * len(hoscflist)) rightb = int(pickwindow[1] / self.getCut()[1] * len(hoscflist))
# threshold = folm * max(hoscflist[leftb : rightb]) # combination of local maximum and threshold
### TEST TEST
threshold = folm * (max(hoscflist[leftb: rightb]) - min(hoscflist[leftb: rightb])) + min( threshold = folm * (max(hoscflist[leftb: rightb]) - min(hoscflist[leftb: rightb])) + min(
hoscflist[leftb: rightb]) # combination of local maximum and threshold hoscflist[leftb: rightb]) # combination of local maximum and threshold
### TEST TEST
m = leftb m = leftb
@ -411,7 +466,6 @@ class SeismicShot(object):
raise ValueError("Distance is NaN for traceID %s" % traceID) raise ValueError("Distance is NaN for traceID %s" % traceID)
return dist return dist
# return abs(float(self.getSrcLoc(traceID))-float(self.getRecLoc(traceID)))
def getRecLoc(self, traceID): ########## input FILENAME ########## def getRecLoc(self, traceID): ########## input FILENAME ##########
''' '''
@ -432,9 +486,7 @@ class SeismicShot(object):
z = coordlist[i].split()[3] z = coordlist[i].split()[3]
return float(x), float(y), float(z) return float(x), float(y), float(z)
# print "WARNING: traceID %s not found" % traceID
raise ValueError("traceID %s not found" % traceID) raise ValueError("traceID %s not found" % traceID)
# return float(self.getSingleStream(traceID)[0].stats.seg2['RECEIVER_LOCATION'])
def getSrcLoc(self): ########## input FILENAME ########## def getSrcLoc(self): ########## input FILENAME ##########
''' '''
@ -477,26 +529,6 @@ class SeismicShot(object):
if len(traceID_list) > 0: if len(traceID_list) > 0:
return traceID_list return traceID_list
# def setManualPicks(self, traceID, picklist): ########## picklist momentan nicht allgemein, nur testweise benutzt ##########
# '''
# Sets the manual picks for a receiver with the ID == traceID for comparison.
# :param: traceID
# :type: int
# :param: picklist, list containing the manual picks (mostlikely, earliest, latest).
# :type: list
# '''
# picks = picklist[traceID - 1].split()
# mostlikely = float(picks[1])
# earliest = float(picks[2])
# latest = float(picks[3])
# if not self.manualpicks.has_key(traceID):
# self.manualpicks[traceID] = (mostlikely, earliest, latest)
# else:
# raise KeyError('MANUAL pick to be set more than once for traceID %s' % traceID)
def setManualPicksFromFile(self, directory='picks'): def setManualPicksFromFile(self, directory='picks'):
''' '''
Read manual picks from *.pck file. Read manual picks from *.pck file.
@ -565,7 +597,6 @@ class SeismicShot(object):
:param: (tnoise, tgap, tsignal), as used in pylot SNR :param: (tnoise, tgap, tsignal), as used in pylot SNR
''' '''
from pylot.core.pick.utils import getSNR from pylot.core.pick.utils import getSNR
tgap = self.getTgap() tgap = self.getTgap()
@ -591,7 +622,7 @@ class SeismicShot(object):
return distancearray return distancearray
def plot2dttc(self, ax=None): ########## 2D ########## def plot2dttc(self, ax=None):
''' '''
Function to plot the traveltime curve for automated picks of a shot. 2d only! ATM: X DIRECTION!! Function to plot the traveltime curve for automated picks of a shot. 2d only! ATM: X DIRECTION!!
''' '''
@ -617,7 +648,7 @@ class SeismicShot(object):
ax.text(0.5, 0.9, 'shot: %s' % self.getShotnumber(), transform=ax.transAxes ax.text(0.5, 0.9, 'shot: %s' % self.getShotnumber(), transform=ax.transAxes
, horizontalalignment='center') , horizontalalignment='center')
def plotmanual2dttc(self, ax=None): ########## 2D ########## def plotmanual2dttc(self, ax=None):
''' '''
Function to plot the traveltime curve for manual picks of a shot. 2D only! Function to plot the traveltime curve for manual picks of a shot. 2D only!
''' '''
@ -643,24 +674,6 @@ class SeismicShot(object):
y.append(point[1]) y.append(point[1])
ax.plot(x, y, 'b', label="Manual Picks") ax.plot(x, y, 'b', label="Manual Picks")
# def plotpickwindow(self): ########## 2D ##########
# '''
# Plots the pickwindow of a shot for the 2nd iteration step of picking. 2D only!
# '''
# import matplotlib.pyplot as plt
# plt.interactive('True')
# pickwindowarray_lowerb = []
# pickwindowarray_upperb = []
# i = 1
# for traceID in self.pwindow.keys():
# pickwindowarray_lowerb.append(self.pwindow[traceID][0])
# pickwindowarray_upperb.append(self.pwindow[traceID][1])
# i += 1
# plt.plot(self.getDistArray4ttcPlot(), pickwindowarray_lowerb, ':k')
# plt.plot(self.getDistArray4ttcPlot(), pickwindowarray_upperb, ':k')
def plotTrace(self, traceID, plotSNR=True, lw=1): def plotTrace(self, traceID, plotSNR=True, lw=1):
fig = plt.figure() fig = plt.figure()
ax = fig.add_subplot(111) ax = fig.add_subplot(111)
@ -679,7 +692,7 @@ class SeismicShot(object):
ax.legend() ax.legend()
ax.text(0.05, 0.9, 'SNR: %s' % snr, transform=ax.transAxes) ax.text(0.05, 0.9, 'SNR: %s' % snr, transform=ax.transAxes)
def plot_traces(self, traceID): ########## 2D, muss noch mehr verbessert werden ########## def plot_traces(self, traceID):
from matplotlib.widgets import Button from matplotlib.widgets import Button
def onclick(event): def onclick(event):

72
pylot/core/active/surveyUtils.py
View File

@ -16,26 +16,13 @@ def readParameters(parfile, parameter):
return value return value
def setArtificialPick(shot_dict, traceID, pick):
"""
:param shot_dict:
:param traceID:
:param pick:
:return:
"""
for shot in shot_dict.values():
shot.setPick(traceID, pick)
shot.setPickwindow(traceID, shot.getCut())
def fitSNR4dist(shot_dict, shiftdist=30, shiftSNR=100): def fitSNR4dist(shot_dict, shiftdist=30, shiftSNR=100):
""" """
Approach to fit the decreasing SNR with wave travel distance.
:param shot_dict: :param shot_dict: dictionary containing Seismicshot objects (e.g. survey.getShotDict())
:param shiftdist: :param shiftdist: shift compensating curve by a certain distance to the left
:param shiftSNR: :param shiftSNR: shift compensating curve by a certain SNR value to the bottom
:return: :return:
""" """
import numpy as np import numpy as np
@ -87,6 +74,8 @@ def plotFittedSNR(dists, snrthresholds, snrs, snrBestFit):
def setDynamicFittedSNR(shot_dict, shiftdist=30, shiftSNR=100, p1=0.004, p2=-0.0007): def setDynamicFittedSNR(shot_dict, shiftdist=30, shiftSNR=100, p1=0.004, p2=-0.0007):
""" """
Set SNR values for a dictionary containing Seismicshots (e.g. survey.getShotDict())
by parameters calulated from fitSNR4dist.
:param shot_dict: :param shot_dict:
:type shot_dict: dict :type shot_dict: dict
@ -119,6 +108,7 @@ def setDynamicFittedSNR(shot_dict, shiftdist=30, shiftSNR=100, p1=0.004, p2=-0.0
def setConstantSNR(shot_dict, snrthreshold=2.5): def setConstantSNR(shot_dict, snrthreshold=2.5):
""" """
Set a constant SNR value to all Seismicshots in a dictionary (e.g. survey.getShotDict()).
:param shot_dict: :param shot_dict:
:param snrthreshold: :param snrthreshold:
@ -158,42 +148,18 @@ def findTracesInRanges(shot_dict, distancebin, pickbin):
def cleanUp(survey): def cleanUp(survey):
""" """
Cleans up a Survey object by removing frontend information
:param survey: which can not be saved in the pickle format.
:return:
""" """
for shot in survey.data.values(): for shot in survey.data.values():
shot.traces4plot = {} shot.traces4plot = {}
def plotScatterStats4Shots(survey, variable):
# def plotScatterStats(survey, key, ax = None):
# import matplotlib.pyplot as plt
# x = []; y = []; value = []
# stats = survey.getStats()
# for shotnumber in stats.keys():
# if type(value) == list:
# value.append(stats[shotnumber][key][0])
# else:
# value.append(stats[shotnumber][key])
# x.append(survey.data[shotnumber].getSrcLoc()[0])
# y.append(survey.data[shotnumber].getSrcLoc()[1])
# if ax == None:
# fig = plt.figure()
# ax = fig.add_subplot(111)
# sc = ax.scatter(x, y, s = value, c = value)
# plt.xlabel('X')
# plt.ylabel('Y')
# cbar = plt.colorbar(sc)
# cbar.set_label(key)
def plotScatterStats4Shots(survey, key):
""" """
Statistics, scatter plot. Statistics, scatter plot.
key can be 'mean SNR', 'median SNR', 'mean SPE', 'median SPE', or 'picked traces'
:param survey: :param survey:
:param key: :param variable: can be 'mean SNR', 'median SNR', 'mean SPE', 'median SPE', or 'picked traces'
:return: :return:
""" """
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
@ -225,7 +191,7 @@ def plotScatterStats4Shots(survey, key):
for shot in statsShot.keys(): for shot in statsShot.keys():
x.append(statsShot[shot]['x']) x.append(statsShot[shot]['x'])
y.append(statsShot[shot]['y']) y.append(statsShot[shot]['y'])
value.append(statsShot[shot][key]) value.append(statsShot[shot][variable])
fig = plt.figure() fig = plt.figure()
ax = fig.add_subplot(111) ax = fig.add_subplot(111)
@ -239,19 +205,19 @@ def plotScatterStats4Shots(survey, key):
plt.xlabel('X') plt.xlabel('X')
plt.ylabel('Y') plt.ylabel('Y')
cbar = plt.colorbar(sc) cbar = plt.colorbar(sc)
cbar.set_label(key) cbar.set_label(variable)
for shot in statsShot.keys(): for shot in statsShot.keys():
ax.annotate(' %s' % shot.getShotnumber(), xy=(shot.getSrcLoc()[0], shot.getSrcLoc()[1]), ax.annotate(' %s' % shot.getShotnumber(), xy=(shot.getSrcLoc()[0], shot.getSrcLoc()[1]),
fontsize='x-small', color='k') fontsize='x-small', color='k')
def plotScatterStats4Receivers(survey, key): def plotScatterStats4Receivers(survey, variable):
""" """
Statistics, scatter plot. Statistics, scatter plot.
key can be 'mean SNR', 'median SNR', 'mean SPE', 'median SPE', or 'picked traces'
:param survey: :param survey:
:param key: :param variable: can be 'mean SNR', 'median SNR', 'mean SPE', 'median SPE', or 'picked traces'
:return: :return:
""" """
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
@ -283,7 +249,7 @@ def plotScatterStats4Receivers(survey, key):
for traceID in statsRec.keys(): for traceID in statsRec.keys():
x.append(statsRec[traceID]['x']) x.append(statsRec[traceID]['x'])
y.append(statsRec[traceID]['y']) y.append(statsRec[traceID]['y'])
value.append(statsRec[traceID][key]) value.append(statsRec[traceID][variable])
fig = plt.figure() fig = plt.figure()
ax = fig.add_subplot(111) ax = fig.add_subplot(111)
@ -297,7 +263,7 @@ def plotScatterStats4Receivers(survey, key):
plt.xlabel('X') plt.xlabel('X')
plt.ylabel('Y') plt.ylabel('Y')
cbar = plt.colorbar(sc) cbar = plt.colorbar(sc)
cbar.set_label(key) cbar.set_label(variable)
shot = survey.data.values()[0] shot = survey.data.values()[0]
for traceID in shot.getTraceIDlist(): for traceID in shot.getTraceIDlist():

19
pylot/core/util/utils.py
View File

@ -9,7 +9,25 @@ import re
import subprocess import subprocess
from obspy.core import UTCDateTime from obspy.core import UTCDateTime
def _pickle_method(m):
if m.im_self is None:
return getattr, (m.im_class, m.im_func.func_name)
else:
return getattr, (m.im_self, m.im_func.func_name)
def worker(func, input, cores = 'max', async = False):
return result
import multiprocessing
if cores == 'max':
cores = multiprocessing.cpu_count()
pool = multiprocessing.Pool(cores)
if async == True:
result = pool.map_async(func, input)
else:
result = pool.map(func, input)
pool.close()
def demeanTrace(trace, window): def demeanTrace(trace, window):
""" """
returns the DATA where each trace is demean by the average value within returns the DATA where each trace is demean by the average value within
@ -231,7 +249,6 @@ def runProgram(cmd, parameter=None):
output = subprocess.check_output('{} | tee /dev/stderr'.format(cmd), output = subprocess.check_output('{} | tee /dev/stderr'.format(cmd),
shell=True) shell=True)
if __name__ == "__main__": if __name__ == "__main__":
import doctest import doctest