Merge branch 'develop' of ariadne.geophysik.ruhr-uni-bochum.de:/data/git/pylot into develop

This commit is contained in:
Sebastian Wehling-Benatelli 2015-10-20 12:02:01 +02:00
commit 0064ff1889
10 changed files with 656 additions and 980 deletions

View File

@ -19,6 +19,7 @@ class Survey(object):
self.setParametersForShots()
self._removeAllEmptyTraces()
self._updateShots()
self.setArtificialPick(0, 0) # artificial pick at source origin
def _generateSurvey(self):
from obspy.core import read
@ -34,17 +35,23 @@ class Survey(object):
shot_dict[shotnumber] = seismicshot.SeismicShot(obsfile)
shot_dict[shotnumber].setParameters('shotnumber', shotnumber)
self.setArtificialPick(0, 0) # artificial pick at source origin
self.data = shot_dict
print ("Generated Survey object for %d shots" % len(shotlist))
print ("Total number of traces: %d \n" %self.countAllTraces())
def setArtificialPick(self, traceID, pick):
'''
Sets an artificial pick for a traceID of all shots in the survey object.
(This can be used to create a pick with t = 0 at the source origin)
'''
for shot in self.data.values():
shot.setPick(traceID, pick)
def setParametersForShots(self, cutwindow = (0, 0.2), tmovwind = 0.3, tsignal = 0.03, tgap = 0.0007):
if (cutwindow == (0, 0.2) and tmovwind == 0.3 and
tsignal == 0.03 and tgap == 0.0007):
print ("Warning: Standard values used for "
"setParamters. This may not be clever.")
"setParamters. This might not be clever.")
# CHANGE this later. Parameters only needed for survey, not for each shot.
for shot in self.data.values():
shot.setCut(cutwindow)
@ -121,14 +128,13 @@ class Survey(object):
shot.setPickwindow(traceID, pickwin_used)
shot.pickTraces(traceID, windowsize, folm, HosAic) # picker
# ++ TEST: set and check SNR before adding to distance bin ############################
shot.setSNR(traceID)
#if shot.getSNR(traceID)[0] < snrthreshold:
if shot.getSNR(traceID)[0] < shot.getSNRthreshold(traceID):
shot.removePick(traceID)
# -- TEST: set and check SNR before adding to distance bin ############################
if shot.getPick(traceID) is not None:
# set epp and lpp if SNR > 1 (else earllatepicker cant set values)
if shot.getSNR(traceID)[0] > 1:
shot.setEarllatepick(traceID)
tpicksum += (datetime.now() - tstartpick); tpick = tpicksum/count
@ -138,6 +144,22 @@ class Survey(object):
self._update_progress(shot.getShotname(), tend, progress)
print('\npickAllShots: Finished\n')
def recover(self):
'''
Recovers all (accidently) removed picks. Still regards SNR threshold.
'''
print('Recovering survey...')
numpicks = 0
for shot in self.data.values():
for traceID in shot.getTraceIDlist():
if shot.getFlag(traceID) == 0:
shot.setFlag(traceID, 1)
if shot.getSNR(traceID)[0] < shot.getSNRthreshold(traceID):
shot.removePick(traceID)
else:
numpicks += 1
print('Recovered %d picks'%numpicks)
def setArtificialPick(self, traceID, pick):
for shot in self.data.values():
shot.setPick(traceID, pick)
@ -195,7 +217,7 @@ class Survey(object):
for traceID in shot.getTraceIDlist():
snrlist.append(shot.getSNR(traceID)[0])
dist.append(shot.getDistance(traceID))
if shot.getPick(traceID) is not None:
if shot.getFlag(traceID) is not 0:
pickedTraces += 1
info_dict[shot.getShotnumber()] = {'numtraces': numtraces,
'picked traces': [pickedTraces,
@ -236,7 +258,7 @@ class Survey(object):
traceIDlist.sort()
ttfile.writelines(str(self.countPickedTraces(shot)) + '\n')
for traceID in traceIDlist:
if shot.getPick(traceID) is not None:
if shot.getFlag(traceID) is not 0:
pick = shot.getPick(traceID) * fmtomo_factor
delta = shot.getPickError(traceID) * fmtomo_factor
(x, y, z) = shot.getRecLoc(traceID)
@ -253,14 +275,87 @@ class Survey(object):
def countPickedTraces(self, shot):
count = 0
for traceID in shot.getTraceIDlist():
if shot.getPick(traceID) is not None:
if shot.getFlag(traceID) is not 0:
count += 1
return count
def plotAllPicks(self, plotDeleted = False):
def countAllPickedTraces(self):
count = 0
for shot in self.data.values():
for traceID in shot.getTraceIDlist():
if shot.getFlag(traceID) is not 0:
count += 1
return count
def plotAllShots(self, rows = 3, columns = 4):
'''
Plots all picks over the distance between source and receiver. Returns (ax, region)
Plots all shots as Matrices with the color corresponding to the traveltime for each receiver.
IMPORTANT NOTE: Topography (z - coordinate) is not considered in the diagrams!
'''
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
plt.interactive(True)
fig = plt.figure()
ax = fig.add_subplot(111)
figPerSubplot = columns * rows
index = 1
#shotnames = []
#shotnumbers = []
# for shot in self.data.values():
# shotnames.append(shot.getShotname())
# shotnumbers.append(shot.getShotnumber())
# shotnumbers = [shotnumbers for (shotnumbers, shotnames) in sorted(zip(shotnumbers, shotnames))]
for shotnumber in self.getShotlist():
if index <= figPerSubplot:
#ax = fig.add_subplot(3,3,i, projection = '3d', title = 'shot:'
#+str(shot_dict[shotnumber].getShotnumber()), xlabel = 'X', ylabel = 'Y', zlabel = 'traveltime')
#shot_dict[shotnumber].plot3dttc(ax = ax, plotpicks = True)
ax = fig.add_subplot(3, 4, index)
self.getShot(shotnumber).matshow(ax = ax, colorbar = False, annotations = True)
index += 1
if index > figPerSubplot:
fig.subplots_adjust(left = 0, bottom = 0, right = 1, top = 1, wspace = 0, hspace = 0)
fig = plt.figure()
index = 1
fig.subplots_adjust(left = 0, bottom = 0, right = 1, top = 1, wspace = 0, hspace = 0)
def plotAllPicks(self, plotRemoved = False, colorByVal = 'log10SNR', ax = None, refreshPlot = False):
'''
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)
:param: plotRemoved, if True plots traces that were picked but removed from the survey (flag = 0)
:type: logical
:param: colorByVal, can be "log10SNR", "pickerror", or "spe"
:type: str
Examples:
regions.chooseRectangles():
- lets the user choose several rectangular regions in the plot
regions.plotTracesInRegions():
- creates plots (shot.plot_traces) for all traces in the active regions (i.e. chosen by e.g. chooseRectangles)
regions.setActiveRegionsForDeletion():
- highlights all shots in a the active regions for deletion
regions.deleteMarkedPicks():
- deletes the picks (pick flag set to 0) for all shots set for deletion
regions.deselectSelection(number):
- deselects the region of number = number
'''
import matplotlib.pyplot as plt
import math
plt.interactive(True)
@ -268,30 +363,40 @@ class Survey(object):
dist = []
pick = []
snrloglist = []
snrlog = []
pickerror = []
spe = []
for shot in self.data.values():
for traceID in shot.getTraceIDlist():
if plotDeleted == False:
if shot.getPick(traceID) is not None:
if plotRemoved == False:
if shot.getFlag(traceID) is not 0 or plotRemoved == True:
dist.append(shot.getDistance(traceID))
pick.append(shot.getPick(traceID))
snrloglist.append(math.log10(shot.getSNR(traceID)[0]))
elif plotDeleted == True:
dist.append(shot.getDistance(traceID))
pick.append(shot.getPick(traceID))
snrloglist.append(math.log10(shot.getSNR(traceID)[0]))
snrlog.append(math.log10(shot.getSNR(traceID)[0]))
pickerror.append(shot.getPickError(traceID))
spe.append(shot.getSymmetricPickError(traceID))
ax = self.createPlot(dist, pick, snrloglist, label = 'log10(SNR)')
region = regions(ax, self.data)
color = {'log10SNR': snrlog,
'pickerror': pickerror,
'spe': spe}
if refreshPlot is False:
ax = self.createPlot(dist, pick, color[colorByVal], label = '%s'%colorByVal)
region = regions(ax, self)
ax.legend()
return ax, region
elif refreshPlot is True:
ax = self.createPlot(dist, pick, color[colorByVal], label = '%s'%colorByVal, ax = ax)
ax.legend()
return ax
def createPlot(self, dist, pick, inkByVal, label):
def createPlot(self, dist, pick, inkByVal, label, ax = None):
import matplotlib.pyplot as plt
plt.interactive(True)
cm = plt.cm.jet
if ax is None:
print('Generating new plot...')
fig = plt.figure()
ax = fig.add_subplot(111)
fig = ax.scatter(dist, pick, cmap = cm, c = inkByVal, s = 5, edgecolors = 'none', label = label)
@ -300,6 +405,8 @@ class Survey(object):
plt.title('Plot of all Picks')
plt.xlabel('Distance [m]')
plt.ylabel('Time [s]')
else:
ax.scatter(dist, pick, cmap = cm, c = inkByVal, s = 5, edgecolors = 'none', label = label)
return ax

View File

@ -128,7 +128,12 @@ def rays2VTK(fnin, fdirout = './vtk_files/', nthPoint = 50):
raynumber += 1
firstline = infile.readline()
if firstline == '': break # break at EOF
raynumber = int(firstline.split()[0])
shotnumber = int(firstline.split()[1])
rayValid = int(firstline.split()[4]) # is zero if the ray is invalid
if rayValid == 0:
print('Invalid ray number %d for shot number %d'%(raynumber, shotnumber))
continue
nRayPoints = int(infile.readline().split()[0])
if not shotnumber in rays.keys():
rays[shotnumber] = {}

View File

@ -1,107 +0,0 @@
# -*- coding: utf-8 -*-
import sys
from obspy import read
from obspy import Stream
from obspy import Trace
from datetime import datetime
import numpy as np
from pylot.core.active import surveyUtils
from pylot.core.active import seismicshot
import activeSeismoPick
reload(seismicshot)
reload(surveyUtils)
#####################################################################################
# parameter definitions:#############################################################
#traceslist = list(np.arange(1, 49, 1)) # traces (1-48)
#shotlist = list(np.arange(302, 352, 1)) # shot-files(.dat) (Paffrath: 302-352) (Hauburg: 353-401) (arange+1!)
cutwindow = (0, 0.2) # cut out a part of the trace [seconds]
tmovwind = 0.3 # size of the moving window
windowsize = (5, 0) # windowsize for AIC picker (indices around HOS picks!!)
pickwindow = cutwindow # for local max and threshold picker: fraction of the seismogram used (0...1) TO BE DONE: depends on cutwindow!!!!
folm = 0.6
rockeskyll = False
if rockeskyll == True:
receiverfile = "Geophone_interpoliert_rockes"
sourcefile = "Schusspunkte_rockes"
obsdir = "/rscratch/minos22/marcel/flachseismik/rockeskyll_200615_270615/"
filename = 'survey_rockes.pickle'
else:
receiverfile = "Geophone_interpoliert_GZB"
sourcefile = "Schusspunkte_GZB"
obsdir = "/rscratch/minos22/marcel/flachseismik/GZB_26_06_15_01/"
filename = 'survey_GZB.pickle'
# SNR
tsignal = 0.03
tgap = 0.0007
snrthreshold = 1
######################################################################################
vmin = 333
vmax = 5500
distBinSize = 2
###########################################
############# Settings: ###################
thresholdpick=True
localmaxpick=False
if thresholdpick == True: pickmethod = "threshold"
if localmaxpick == True: pickmethod = "localmax"
HosAic = 'hos' # can be 'hos' or 'aic'
###########################################
starttime = datetime.now()
print '\n--------------------Starting Script at %s -------------------\n' %starttime.time()
if thresholdpick == True: print 'Using treshold pickmethod!\n'
elif localmaxpick == True: print 'Using local maximum pick method!\n'
if HosAic == 'aic': print 'picking with AIC\n'
if HosAic == 'hos': print 'picking with HOS\n'
survey = activeSeismoPick.Survey(obsdir, sourcefile, receiverfile, True)
surveyUtils.setFittedSNR(survey.getShotDict())
print '\nDone after %s seconds!\n------------------------------------------------------------------------------\n' % (datetime.now() - starttime).seconds
count = 0; tpicksum = starttime - starttime
for shot in survey.data.values():
tstartpick = datetime.now(); count += 1
for traceID in shot.getTraceIDlist():
distance = shot.getDistance(traceID) # receive distance
pickwin_used = pickwindow # use pickwindow set in the parameter section
# for higher distances use a linear vmin/vmax to cut out late/early regions with high noise
if distance > 5.:
pwleft = distance/vmax ################## TEST
pwright = distance/vmin
if pwright > cutwindow[1]:
pwright = cutwindow[1]
pickwin_used = (pwleft, pwright)
shot.setPickwindow(traceID, pickwin_used)
shot.pickTraces(traceID, windowsize, folm, HosAic) # picker
#shot.setManualPicks(traceID, picklist) # set manual picks if given (yet used on 2D only)
# ++ TEST: set and check SNR before adding to distance bin ############################
shot.setSNR(traceID)
#if shot.getSNR(traceID)[0] < snrthreshold:
if shot.getSNR(traceID)[0] < shot.getSNRthreshold(traceID):
shot.removePick(traceID)
# -- TEST: set and check SNR before adding to distance bin ############################
if shot.getPick(traceID) is not None:
shot.setEarllatepick(traceID)
tpicksum += (datetime.now() - tstartpick); tpick = tpicksum/count
tremain = (tpick * (len(survey.getShotDict()) - count))
tend = datetime.now() + tremain
print 'shot: %s, est. time to be finished is %s:%s:%s' % (shot.getShotname(), tend.hour, tend.minute, tend.second)
survey.saveSurvey(filename)
print '\n--- Finished script ---'
print 'Elapsed time:', datetime.now()-starttime

View File

@ -491,7 +491,11 @@ class SeisArray(object):
plt.legend()
if annotations == True:
for traceID in self.getReceiverCoordinates().keys():
plt.annotate(str(traceID), xy = (self._getXreceiver(traceID), self._getYreceiver(traceID)), fontsize = 'x-small')
plt.annotate((' ' + str(traceID)), xy = (self._getXreceiver(traceID), self._getYreceiver(traceID)), fontsize = 'x-small', color = 'k')
for shotnumber in self.getSourceLocations().keys():
plt.annotate((' ' + str(shotnumber)), xy = (self._getXshot(shotnumber), self._getYshot(shotnumber)), fontsize = 'x-small', color = 'b')
def plotArray3D(self, ax = None):
import matplotlib.pyplot as plt

View File

@ -1,666 +0,0 @@
# -*- coding: utf-8 -*-
import sys
import numpy as np
from scipy.interpolate import griddata
class SeisArray(object):
'''
Can be used to interpolate missing values of a receiver grid, if only support points were measured.
Input file should contain in each line: ('traceID' 'receiverLineID' 'number of the geophone on recLine' 'X' 'Y' 'Z')
Can be used to generate a velocity grid file (vgrids.in) for FMTOMO with a topography adapting gradient.
Can be used to generate an interface file for FMTOMO (right now only interface.z used by grid3dg) for the topography.
Supports vtk output for sources and receivers.
Note: Source and Receiver files for FMTOMO will be generated by the Survey object (because traveltimes will be added directly).
'''
def __init__(self, recfile):
self._receiverlines = {}
self._receiverCoords = {}
self._measuredReceivers = {}
self._measuredTopo = {}
self._sourceLocs = {}
self._geophoneNumbers = {}
self._receiverlist = open(recfile, 'r').readlines()
self._generateReceiverlines()
self._setReceiverCoords()
self._setGeophoneNumbers()
def _generateReceiverlines(self):
'''
Connects the traceIDs to the lineIDs
for each receiverline in a dictionary.
'''
for receiver in self._receiverlist:
traceID = int(receiver.split()[0])
lineID = int(receiver.split()[1])
if not lineID in self._receiverlines.keys():
self._receiverlines[lineID] = []
self._receiverlines[lineID].append(traceID)
def _setReceiverCoords(self):
'''
Fills the three x, y, z dictionaries with measured coordinates
'''
for line in self._getReceiverlist():
traceID = int(line.split()[0])
x = float(line.split()[3])
y = float(line.split()[4])
z = float(line.split()[5])
self._receiverCoords[traceID] = (x, y, z)
self._measuredReceivers[traceID] = (x, y, z)
def _setGeophoneNumbers(self):
for line in self._getReceiverlist():
traceID = int(line.split()[0])
gphoneNum = float(line.split()[2])
self._geophoneNumbers[traceID] = gphoneNum
def _getReceiverlines(self):
return self._receiverlines
def _getReceiverlist(self):
return self._receiverlist
def getReceiverCoordinates(self):
return self._receiverCoords
def _getXreceiver(self, traceID):
return self._receiverCoords[traceID][0]
def _getYreceiver(self, traceID):
return self._receiverCoords[traceID][1]
def _getZreceiver(self, traceID):
return self._receiverCoords[traceID][2]
def _getXshot(self, shotnumber):
return self._sourceLocs[shotnumber][0]
def _getYshot(self, shotnumber):
return self._sourceLocs[shotnumber][1]
def _getZshot(self, shotnumber):
return self._sourceLocs[shotnumber][2]
def _getReceiverValue(self, traceID, coordinate):
setCoordinate = {'X': self._getXreceiver,
'Y': self._getYreceiver,
'Z': self._getZreceiver}
return setCoordinate[coordinate](traceID)
def _getGeophoneNumber(self, traceID):
return self._geophoneNumbers[traceID]
def getMeasuredReceivers(self):
return self._measuredReceivers
def getMeasuredTopo(self):
return self._measuredTopo
def getSourceLocations(self):
return self._sourceLocs
def _setXvalue(self, traceID, value):
self._checkKey(traceID)
self._receiverCoords[traceID][0] = value
def _setYvalue(self, traceID, value):
self._checkKey(traceID)
self._receiverCoords[traceID][1] = value
def _setZvalue(self, traceID, value):
self._checkKey(traceID)
self._receiverCoords[traceID][2] = value
def _setValue(self, traceID, coordinate, value):
setCoordinate = {'X': self._setXvalue,
'Y': self._setYvalue,
'Z': self._setZvalue}
setCoordinate[coordinate](traceID, value)
def _checkKey(self, traceID):
if not traceID in self._receiverCoords.keys():
self._receiverCoords[traceID] = [None, None, None]
def _checkTraceIDdirection(self, traceID1, traceID2):
if traceID2 > traceID1:
direction = +1
return direction
if traceID2 < traceID1:
direction = -1
return direction
print "Error: Same Value for traceID1 = %s and traceID2 = %s" %(traceID1, traceID2)
def _checkCoordDirection(self, traceID1, traceID2, coordinate):
'''
Checks whether the interpolation direction is positive or negative
'''
if self._getReceiverValue(traceID1, coordinate) < self._getReceiverValue(traceID2, coordinate):
direction = +1
return direction
if self._getReceiverValue(traceID1, coordinate) > self._getReceiverValue(traceID2, coordinate):
direction = -1
return direction
print "Error: Same Value for traceID1 = %s and traceID2 = %s" %(traceID1, traceID2)
def _interpolateMeanDistances(self, traceID1, traceID2, coordinate):
'''
Returns the mean distance between two traceID's depending on the number of geophones in between
'''
num_spaces = abs(self._getGeophoneNumber(traceID1) - self._getGeophoneNumber(traceID2))
mean_distance = abs(self._getReceiverValue(traceID1, coordinate) - self._getReceiverValue(traceID2, coordinate))/num_spaces
return mean_distance
def interpolateValues(self, coordinate):
'''
Interpolates and sets all values (linear) for coordinate = 'X', 'Y' or 'Z'
'''
for lineID in self._getReceiverlines().keys():
number_measured = len(self._getReceiverlines()[lineID])
for index, traceID1 in enumerate(self._getReceiverlines()[lineID]):
if index + 1 < number_measured:
traceID2 = self._getReceiverlines()[lineID][index + 1]
traceID_dir = self._checkTraceIDdirection(traceID1, traceID2)
traceID_interp = traceID1 + traceID_dir
coord_dir = self._checkCoordDirection(traceID1, traceID2, coordinate)
mean_distance = self._interpolateMeanDistances(traceID1, traceID2, coordinate)
while (traceID_dir * traceID_interp) < (traceID_dir * traceID2):
self._setValue(traceID_interp, coordinate,
(self._getReceiverValue(traceID_interp - traceID_dir, coordinate)
+ coord_dir * (mean_distance)))
traceID_interp += traceID_dir
def addMeasuredTopographyPoints(self, filename):
'''
Use more measured points for a higher precision of height interpolation.
Input file should contain in each line: ('point ID' 'X' 'Y' 'Z')
'''
topolist = open(filename, 'r').readlines()
for line in topolist:
line = line.split()
pointID = int(line[0])
x = float(line[1])
y = float(line[2])
z = float(line[3])
self._measuredTopo[pointID] = (x, y, z)
def addSourceLocations(self, filename):
'''
Use source locations for a higher precision of height interpolation.
Input file should contain in each line: ('point ID' 'X' 'Y' 'Z')
Source locations must be added before they can be written to vtk files.
'''
topolist = open(filename, 'r').readlines()
for line in topolist:
line = line.split()
pointID = int(line[0])
x = float(line[1])
y = float(line[2])
z = float(line[3])
self._sourceLocs[pointID] = (x, y, z)
def interpZcoords4rec(self, method = 'linear'):
'''
Interpolates z values for all receivers.
'''
measured_x, measured_y, measured_z = self.getAllMeasuredPointsLists()
for traceID in self.getReceiverCoordinates().keys():
if type(self.getReceiverCoordinates()[traceID]) is not tuple:
z = griddata((measured_x, measured_y), measured_z, (self._getXreceiver(traceID), self._getYreceiver(traceID)), method = method)
self._setZvalue(traceID, float(z))
def _getAngle(self, distance):
'''
Function returns the angle on a Sphere of the radius R = 6371 [km] for a distance [km].
'''
PI = np.pi
R = 6371.
angle = distance * 180 / (PI * R)
return angle
def _getDistance(self, angle):
'''
Function returns the distance [km] on a Sphere of the radius R = 6371 [km] for an angle.
'''
PI = np.pi
R = 6371.
distance = angle / 180 * (PI * R)
return distance
def getMeasuredReceiverLists(self):
'''
Returns a list of all measured receivers known to SeisArray.
'''
x = []; y = []; z = []
for traceID in self.getMeasuredReceivers().keys():
x.append(self.getMeasuredReceivers()[traceID][0])
y.append(self.getMeasuredReceivers()[traceID][1])
z.append(self.getMeasuredReceivers()[traceID][2])
return x, y, z
def getMeasuredTopoLists(self):
'''
Returns a list of all measured topography points known to the SeisArray.
'''
x = []; y = []; z = []
for pointID in self.getMeasuredTopo().keys():
x.append(self.getMeasuredTopo()[pointID][0])
y.append(self.getMeasuredTopo()[pointID][1])
z.append(self.getMeasuredTopo()[pointID][2])
return x, y, z
def getSourceLocsLists(self):
'''
Returns a list of all measured source locations known to SeisArray.
'''
x = []; y = []; z = []
for pointID in self.getSourceLocations().keys():
x.append(self.getSourceLocations()[pointID][0])
y.append(self.getSourceLocations()[pointID][1])
z.append(self.getSourceLocations()[pointID][2])
return x, y, z
def getAllMeasuredPointsLists(self):
'''
Returns a list of all measured points known to SeisArray.
'''
mtopo_x, mtopo_y, mtopo_z = self.getMeasuredTopoLists()
msource_x, msource_y, msource_z = self.getSourceLocsLists()
mrec_x, mrec_y, mrec_z = self.getMeasuredReceiverLists()
x = mtopo_x + mrec_x + msource_x
y = mtopo_y + mrec_y + msource_y
z = mtopo_z + mrec_z + msource_z
return x, y, z
def getReceiverLists(self):
'''
Returns a list of all receivers (measured and interpolated).
'''
x = []; y =[]; z = []
for traceID in self.getReceiverCoordinates().keys():
x.append(self.getReceiverCoordinates()[traceID][0])
y.append(self.getReceiverCoordinates()[traceID][1])
z.append(self.getReceiverCoordinates()[traceID][2])
return x, y, z
def _interpolateXY4rec(self):
'''
Interpolates the X and Y coordinates for all receivers.
'''
for coordinate in ('X', 'Y'):
self.interpolateValues(coordinate)
def interpolateAll(self):
self._interpolateXY4rec()
self.interpZcoords4rec()
def interpolateTopography(self, nTheta, nPhi, thetaSN, phiWE, method = 'linear', filename = 'interface1.in'):
'''
Interpolate Z values on a regular grid with cushion nodes to use it as FMTOMO topography interface.
Returns a surface in form of a list of points [[x1, y1, z1], [x2, y2, y2], ...] (cartesian).
:param: nTheta, number of points in theta (NS)
type: integer
:param: nPhi, number of points in phi (WE)
type: integer
:param: thetaSN (S, N) extensions of the model in degree
type: tuple
:param: phiWE (W, E) extensions of the model in degree
type: tuple
'''
surface = []
elevation = 0.25 # elevate topography so that no source lies above the surface
if filename is not None:
outfile = open(filename, 'w')
print "Interpolating topography on regular grid with the dimensions:"
print "nTheta = %s, nPhi = %s, thetaSN = %s, phiWE = %s"%(nTheta, nPhi, thetaSN, phiWE)
print "method = %s, filename = %s" %(method, filename)
thetaS, thetaN = thetaSN
phiW, phiE = phiWE
measured_x, measured_y, measured_z = self.getAllMeasuredPointsLists()
# need to determine the delta to add two cushion nodes around the min/max values
thetaDelta = (thetaN - thetaS) / (nTheta - 1)
phiDelta = (phiE - phiW) / (nPhi - 1)
thetaGrid = np.linspace(thetaS - thetaDelta, thetaN + thetaDelta, num = nTheta + 2) # +2 cushion nodes
phiGrid = np.linspace(phiW - phiDelta, phiE + phiDelta, num = nPhi + 2) # +2 cushion nodes
nTotal = len(thetaGrid) * len(phiGrid); count = 0
for theta in thetaGrid:
for phi in phiGrid:
xval = self._getDistance(phi)
yval = self._getDistance(theta)
z = griddata((measured_x, measured_y), measured_z, (xval, yval), method = method)
# in case the point lies outside, nan will be returned. Find nearest:
if np.isnan(z) == True:
z = griddata((measured_x, measured_y), measured_z, (xval, yval), method = 'nearest')
z = float(z)
surface.append((xval, yval, z))
count += 1
progress = float(count) / float(nTotal) * 100
self._update_progress(progress)
if filename is not None:
outfile.writelines('%10s\n'%(z + elevation))
return surface
def generateVgrid(self, nTheta = 80, nPhi = 80, nR = 120,
thetaSN = (-0.2, 1.2), phiWE = (-0.2, 1.2),
Rbt = (-62.0, 6.0), vbot = 5.5, filename = 'vgrids.in',
method = 'linear' ):
'''
Generate a velocity grid for fmtomo regarding topography with a linear gradient starting at the topography with 0.34 [km/s].
:param: nTheta, number of points in theta (NS)
type: integer
:param: nPhi, number of points in phi (WE)
type: integer
:param: nR, number of points in depth
type: integer
:param: thetaSN (S, N) extensions of the model in degree
type: tuple
:param: phiWE (W, E) extensions of the model in degree
type: tuple
:param: Rbt (bot, top) extensions of the model in km
type: tuple
:param: vbot, velocity at the bottom of the model
type: real
'''
def getRad(angle):
PI = np.pi
rad = angle / 180 * PI
return rad
def getZmax(surface):
z = []
for point in surface:
z.append(point[2])
return max(z)
R = 6371
vmin = 0.34
decm = 0.3 # diagonal elements of the covariance matrix (grid3dg's default value is 0.3)
outfile = open(filename, 'w')
thetaS, thetaN = thetaSN
phiW, phiE = phiWE
rbot = Rbt[0] + R
rtop = Rbt[1] + R
# need to determine the delta to add two cushion nodes around the min/max values
thetaDelta = abs(thetaN - thetaS) / (nTheta - 1)
phiDelta = abs(phiE - phiW) / (nPhi - 1)
rDelta = abs(rbot - rtop) / (nR - 1)
# create a regular grid including +2 cushion nodes in every direction
thetaGrid = np.linspace(thetaS - thetaDelta, thetaN + thetaDelta, num = nTheta + 2) # +2 cushion nodes
phiGrid = np.linspace(phiW - phiDelta, phiE + phiDelta, num = nPhi + 2) # +2 cushion nodes
rGrid = np.linspace(rbot - rDelta, rtop + rDelta, num = nR + 2) # +2 cushion nodes
nTotal = len(rGrid) * len(thetaGrid) * len(phiGrid)
print "Total number of grid nodes: %s"%nTotal
# write header for velocity grid file (in RADIANS)
outfile.writelines('%10s %10s \n' %(1, 1))
outfile.writelines('%10s %10s %10s\n' %(nR + 2, nTheta + 2, nPhi + 2))
outfile.writelines('%10s %10s %10s\n' %(rDelta, getRad(thetaDelta), getRad(phiDelta)))
outfile.writelines('%10s %10s %10s\n' %(rbot - rDelta, getRad(thetaS - thetaDelta), getRad(phiW - phiDelta)))
surface = self.interpolateTopography(nTheta, nPhi, thetaSN, phiWE, method = method, filename = None)
zmax = getZmax(surface)
print "\nGenerating velocity grid for FMTOMO. Output filename = %s, interpolation method = %s"%(filename, method)
print "nTheta = %s, nPhi = %s, nR = %s, thetaSN = %s, phiWE = %s, Rbt = %s"%(nTheta, nPhi, nR, thetaSN, phiWE, Rbt)
count = 0
for radius in rGrid:
for theta in thetaGrid:
for phi in phiGrid:
xval = self._getDistance(phi)
yval = self._getDistance(theta)
for point in surface:
if point[0] == xval and point[1] == yval:
z = point[2]
if radius > (R + z + 1):
vel = 0.0
# elif radius > (R + z - 15): ########### TESTING
# vel = (radius - z - R) / (Rbt[0] - rDelta - zmax) * 1.0 + vmin ##########################
else:
vel = (radius - z - R) / (Rbt[0] - rDelta - zmax) * vbot + vmin ##########################
count += 1
outfile.writelines('%10s %10s\n'%(vel, decm))
progress = float(count) / float(nTotal) * 100
self._update_progress(progress)
outfile.close()
def exportAll(self, filename = 'interpolated_receivers.out'):
recfile_out = open(filename, 'w')
count = 0
for traceID in self.getReceiverCoordinates().keys():
count += 1
x, y, z = self.getReceiverCoordinates()[traceID]
recfile_out.writelines('%5s %15s %15s %15s\n' %(traceID, x, y, z))
print "Exported coordinates for %s traces to file > %s" %(count, filename)
recfile_out.close()
def plotArray2D(self, plot_topo = False, highlight_measured = False, annotations = True):
import matplotlib.pyplot as plt
plt.interactive(True)
plt.figure()
xmt, ymt, zmt = self.getMeasuredTopoLists()
xsc, ysc, zsc = self.getSourceLocsLists()
xmr, ymr, zmr = self.getMeasuredReceiverLists()
xrc, yrc, zrc = self.getReceiverLists()
plt.plot(xrc, yrc, 'k.', markersize = 10, label = 'all receivers')
plt.plot(xsc, ysc, 'b*', markersize = 10, label = 'shot locations')
if plot_topo == True:
plt.plot(xmt, ymt, 'b', markersize = 10, label = 'measured topo points')
if highlight_measured == True:
plt.plot(xmr, ymr, 'ro', label = 'measured receivers')
plt.xlabel('X [m]')
plt.ylabel('Y [m]')
plt.legend()
if annotations == True:
for traceID in self.getReceiverCoordinates().keys():
plt.annotate(str(traceID), xy = (self._getXreceiver(traceID), self._getYreceiver(traceID)), fontsize = 'x-small')
def plotArray3D(self, ax = None):
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
plt.interactive(True)
if ax == None:
fig = plt.figure()
ax = plt.axes(projection = '3d')
xmt, ymt, zmt = self.getMeasuredTopoLists()
xmr, ymr, zmr = self.getMeasuredReceiverLists()
xin, yin, zin = self.getReceiverLists()
ax.plot(xmt, ymt, zmt, 'b*', markersize = 10, label = 'measured topo points')
ax.plot(xin, yin, zin, 'k.', markersize = 10, label = 'interpolated receivers')
ax.plot(xmr, ymr, zmr, 'ro', label = 'measured receivers')
ax.set_xlabel('X'); ax.set_ylabel('Y'); ax.set_zlabel('elevation')
ax.legend()
return ax
def plotSurface3D(self, ax = None, step = 0.5, method = 'linear'):
from matplotlib import cm
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
plt.interactive(True)
if ax == None:
fig = plt.figure()
ax = plt.axes(projection = '3d')
xmt, ymt, zmt = self.getMeasuredTopoLists()
xmr, ymr, zmr = self.getMeasuredReceiverLists()
x = xmt + xmr
y = ymt + ymr
z = zmt + zmr
xaxis = np.arange(min(x)+1, max(x), step)
yaxis = np.arange(min(y)+1, max(y), step)
xgrid, ygrid = np.meshgrid(xaxis, yaxis)
zgrid = griddata((x, y), z, (xgrid, ygrid), method = method)
ax.plot_surface(xgrid, ygrid, zgrid, linewidth = 0, cmap = cm.jet, vmin = min(z), vmax = max(z))
ax.set_zlim(-(max(x) - min(x)/2),(max(x) - min(x)/2))
ax.set_aspect('equal')
ax.set_xlabel('X'); ax.set_ylabel('Y'); ax.set_zlabel('elevation')
ax.legend()
return ax
def _update_progress(self, progress):
sys.stdout.write("%d%% done \r" % (progress) )
sys.stdout.flush()
def receivers2VTK(self, filename = 'receivers.vtk'):
'''
Generates vtk files from all receivers of the SeisArray object.
'''
outfile = open(filename, 'w')
traceIDs = []
for traceID in self.getReceiverCoordinates():
traceIDs.append(traceID)
nPoints = len(traceIDs)
# write header
print("Writing header for VTK file...")
outfile.writelines('# vtk DataFile Version 3.1\n')
outfile.writelines('Receivers with traceIDs\n')
outfile.writelines('ASCII\n')
outfile.writelines('DATASET POLYDATA\n')
outfile.writelines('POINTS %15d float\n' %(nPoints))
# write coordinates
print("Writing coordinates to VTK file...")
for traceID in traceIDs:
x = self._getXreceiver(traceID)
y = self._getYreceiver(traceID)
z = self._getZreceiver(traceID)
outfile.writelines('%10f %10f %10f \n' %(x, y, z))
outfile.writelines('VERTICES %15d %15d\n' %(nPoints, 2 * nPoints))
# write indices
print("Writing indices to VTK file...")
for index in range(nPoints):
outfile.writelines('%10d %10d\n' %(1, index))
outfile.writelines('POINT_DATA %15d\n' %(nPoints))
outfile.writelines('SCALARS traceIDs int %d\n' %(1))
outfile.writelines('LOOKUP_TABLE default\n')
# write traceIDs
print("Writing traceIDs to VTK file...")
for traceID in traceIDs:
outfile.writelines('%10d\n' %traceID)
outfile.close()
print("Wrote receiver grid for %d points to file: %s" %(nPoints, filename))
return
def sources2VTK(self, filename = 'sources.vtk'):
'''
Generates vtk-files for all source locations in the SeisArray object.
'''
outfile = open(filename, 'w')
shotnumbers = []
for shotnumber in self.getSourceLocations():
shotnumbers.append(shotnumber)
nPoints = len(shotnumbers)
# write header
print("Writing header for VTK file...")
outfile.writelines('# vtk DataFile Version 3.1\n')
outfile.writelines('Shots with shotnumbers\n')
outfile.writelines('ASCII\n')
outfile.writelines('DATASET POLYDATA\n')
outfile.writelines('POINTS %15d float\n' %(nPoints))
# write coordinates
print("Writing coordinates to VTK file...")
for shotnumber in shotnumbers:
x = self._getXshot(shotnumber)
y = self._getYshot(shotnumber)
z = self._getZshot(shotnumber)
outfile.writelines('%10f %10f %10f \n' %(x, y, z))
outfile.writelines('VERTICES %15d %15d\n' %(nPoints, 2 * nPoints))
# write indices
print("Writing indices to VTK file...")
for index in range(nPoints):
outfile.writelines('%10d %10d\n' %(1, index))
outfile.writelines('POINT_DATA %15d\n' %(nPoints))
outfile.writelines('SCALARS shotnumbers int %d\n' %(1))
outfile.writelines('LOOKUP_TABLE default\n')
# write shotnumber
print("Writing shotnumbers to VTK file...")
for shotnumber in shotnumbers:
outfile.writelines('%10d\n' %shotnumber)
outfile.close()
print("Wrote receiver grid for %d points to file: %s" %(nPoints, filename))
return
def saveSeisArray(self, filename = 'seisArray.pickle'):
import cPickle
outfile = open(filename, 'wb')
cPickle.dump(self, outfile, -1)
print('saved SeisArray to file %s'%(filename))
@staticmethod
def from_pickle(filename):
import cPickle
infile = open(filename, 'rb')
return cPickle.load(infile)

View File

@ -10,6 +10,8 @@ from pylot.core.pick.CharFuns import HOScf
from pylot.core.pick.CharFuns import AICcf
from pylot.core.pick.utils import getSNR
from pylot.core.pick.utils import earllatepicker
import matplotlib.pyplot as plt
plt.interactive('True')
class SeismicShot(object):
'''
@ -27,9 +29,6 @@ class SeismicShot(object):
self.srcCoordlist = None
self.traceIDs = None
self.pick = {}
self.pick_backup = {}
self.earliest = {}
self.latest = {}
self.pickwindow= {}
self.manualpicks= {}
self.snr = {}
@ -132,17 +131,27 @@ class SeismicShot(object):
def getSourcefile(self):
return self.paras['sourcefile']
def getPick(self, traceID):
return self.pick[traceID]
def getPick(self, traceID, returnRemoved = False):
if not self.getFlag(traceID) == 0:
return self.pick[traceID]['mpp']
if returnRemoved == True:
#print('getPick: Returned removed pick for shot %d, traceID %d' %(self.getShotnumber(), traceID))
return self.pick[traceID]['mpp']
def getPick_backup(self, traceID):
return self.pick_backup[traceID]
def getPickIncludeRemoved(self, traceID):
return self.getPick(traceID, returnRemoved = True)
def getEarliest(self, traceID):
return self.earliest[traceID]
return self.pick[traceID]['epp']
def getLatest(self, traceID):
return self.latest[traceID]
return self.pick[traceID]['lpp']
def getSymmetricPickError(self, traceID):
pickerror = self.pick[traceID]['spe']
if np.isnan(pickerror) == True:
print "SPE is NaN for shot %s, traceID %s"%(self.getShotnumber(), traceID)
return pickerror
def getPickError(self, traceID):
pickerror = abs(self.getEarliest(traceID) - self.getLatest(traceID))
@ -217,7 +226,7 @@ class SeismicShot(object):
:type: int
'''
return HOScf(self.getSingleStream(traceID), self.getCut(),
self.getTmovwind(), self.getOrder())
self.getTmovwind(), self.getOrder(), stealthMode = True)
def getAICcf(self, traceID):
'''
@ -240,7 +249,7 @@ class SeismicShot(object):
tr_cf = Trace()
tr_cf.data = self.getHOScf(traceID).getCF()
st_cf += tr_cf
return AICcf(st_cf, self.getCut(), self.getTmovwind())
return AICcf(st_cf, self.getCut(), self.getTmovwind(), stealthMode = True)
def getSingleStream(self, traceID): ########## SEG2 / SEGY ? ##########
'''
@ -305,16 +314,17 @@ class SeismicShot(object):
'aic': aiccftime}
self.setPick(traceID, setHosAic[HosAic])
self.pick_backup[traceID] = setHosAic[HosAic] ### verbessern (vor allem weil ueberschrieben bei 2tem mal picken)
def setEarllatepick(self, traceID, nfac = 1.5):
tgap = self.getTgap()
tsignal = self.getTsignal()
tnoise = self.getPick(traceID) - tgap
tnoise = self.getPickIncludeRemoved(traceID) - tgap
(self.earliest[traceID], self.latest[traceID], tmp) = earllatepicker(self.getSingleStream(traceID),
(self.pick[traceID]['epp'], self.pick[traceID]['lpp'],
self.pick[traceID]['spe']) = earllatepicker(self.getSingleStream(traceID),
nfac, (tnoise, tgap, tsignal),
self.getPick(traceID))
self.getPickIncludeRemoved(traceID),
stealthMode = True)
def threshold(self, hoscf, aiccf, windowsize, pickwindow, folm = 0.6):
'''
@ -464,7 +474,10 @@ class SeismicShot(object):
# raise KeyError('MANUAL pick to be set more than once for traceID %s' % traceID)
def setPick(self, traceID, pick): ########## siehe Kommentar ##########
self.pick[traceID] = pick
if not traceID in self.pick.keys():
self.pick[traceID] = {}
self.pick[traceID]['mpp'] = pick
self.pick[traceID]['flag'] = 1
# ++++++++++++++ Block raus genommen, da Error beim 2ten Mal picken! (Ueberschreiben von erstem Pick!)
# if not self.pick.has_key(traceID):
# self.getPick(traceID) = picks
@ -475,7 +488,14 @@ class SeismicShot(object):
# parlist = open(parfile,'r').readlines()
def removePick(self, traceID):
self.setPick(traceID, None)
self.setFlag(traceID, 0)
def setFlag(self, traceID, flag):
'Set flag = 0 if pick is invalid, else flag = 1'
self.pick[traceID]['flag'] = flag
def getFlag(self, traceID):
return self.pick[traceID]['flag']
def setPickwindow(self, traceID, pickwindow):
self.pickwindow[traceID] = pickwindow
@ -580,7 +600,39 @@ class SeismicShot(object):
# plt.plot(self.getDistArray4ttcPlot(), pickwindowarray_upperb, ':k')
def plot_traces(self, traceID, folm = 0.6): ########## 2D, muss noch mehr verbessert werden ##########
import matplotlib.pyplot as plt
from matplotlib.widgets import Button
def onclick(event):
self.setPick(traceID, event.xdata)
self._drawStream(traceID, refresh = True)
self._drawCFs(traceID, folm, refresh = True)
fig.canvas.mpl_disconnect(self.traces4plot[traceID]['cid'])
plt.draw()
def connectButton(event = None):
cid = fig.canvas.mpl_connect('button_press_event', onclick)
self.traces4plot[traceID]['cid'] = cid
fig = plt.figure()
ax1 = fig.add_subplot(2,1,1)
ax2 = fig.add_subplot(2,1,2, sharex = ax1)
axb = fig.add_axes([0.15, 0.91, 0.05, 0.03])
button = Button(axb, 'repick', color = 'red', hovercolor = 'grey')
button.on_clicked(connectButton)
self.traces4plot = {}
if not traceID in self.traces4plot.keys():
self.traces4plot[traceID] = {'fig': fig,
'ax1': ax1,
'ax2': ax2,
'axb': axb,
'button': button,
'cid': None,}
self._drawStream(traceID)
self._drawCFs(traceID, folm)
def _drawStream(self, traceID, refresh = False):
from pylot.core.util.utils import getGlobalTimes
from pylot.core.util.utils import prepTimeAxis
@ -589,32 +641,47 @@ class SeismicShot(object):
timeaxis = prepTimeAxis(stime, stream[0])
timeaxis -= stime
plt.interactive('True')
ax = self.traces4plot[traceID]['ax1']
hoscf = self.getHOScf(traceID)
aiccf = self.getAICcf(traceID)
if refresh == True:
xlim, ylim = ax.get_xlim(), ax.get_ylim()
ax.clear()
if refresh == True:
ax.set_xlim(xlim)
ax.set_ylim(ylim)
fig = plt.figure()
ax1 = plt.subplot(2,1,1)
plt.title('Shot: %s, traceID: %s, pick: %s' %(self.getShotnumber(), traceID, self.getPick(traceID)))
ax1.plot(timeaxis, stream[0].data, 'k', label = 'trace')
ax1.plot([self.getPick(traceID), self.getPick(traceID)],
ax.set_title('Shot: %s, traceID: %s, pick: %s'
%(self.getShotnumber(), traceID, self.getPick(traceID)))
ax.plot(timeaxis, stream[0].data, 'k', label = 'trace')
ax.plot([self.getPick(traceID), self.getPick(traceID)],
[min(stream[0].data),
max(stream[0].data)],
'r', label = 'mostlikely')
plt.legend()
ax2 = plt.subplot(2,1,2, sharex = ax1)
ax2.plot(hoscf.getTimeArray(), hoscf.getCF(), 'b', label = 'HOS')
ax2.plot(hoscf.getTimeArray(), aiccf.getCF(), 'g', label = 'AIC')
ax2.plot([self.getPick(traceID), self.getPick(traceID)],
ax.legend()
def _drawCFs(self, traceID, folm, refresh = False):
hoscf = self.getHOScf(traceID)
aiccf = self.getAICcf(traceID)
ax = self.traces4plot[traceID]['ax2']
if refresh == True:
xlim, ylim = ax.get_xlim(), ax.get_ylim()
ax.clear()
if refresh == True:
ax.set_xlim(xlim)
ax.set_ylim(ylim)
ax.plot(hoscf.getTimeArray(), hoscf.getCF(), 'b', label = 'HOS')
ax.plot(hoscf.getTimeArray(), aiccf.getCF(), 'g', label = 'AIC')
ax.plot([self.getPick(traceID), self.getPick(traceID)],
[min(np.minimum(hoscf.getCF(), aiccf.getCF())),
max(np.maximum(hoscf.getCF(), aiccf.getCF()))],
'r', label = 'mostlikely')
ax2.plot([0, self.getPick(traceID)],
ax.plot([0, self.getPick(traceID)],
[folm * max(hoscf.getCF()), folm * max(hoscf.getCF())],
'm:', label = 'folm = %s' %folm)
plt.xlabel('Time [s]')
plt.legend()
ax.set_xlabel('Time [s]')
ax.legend()
def plot3dttc(self, step = 0.5, contour = False, plotpicks = False, method = 'linear', ax = None):
'''
@ -632,7 +699,6 @@ class SeismicShot(object):
:param: method (optional), interpolation method; can be 'linear' (default) or 'cubic'
:type: 'string'
'''
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
@ -641,13 +707,13 @@ class SeismicShot(object):
y = []
z = []
for traceID in self.pick.keys():
if self.getPick(traceID) != None:
if self.getFlag(traceID) != 0:
x.append(self.getRecLoc(traceID)[0])
y.append(self.getRecLoc(traceID)[1])
z.append(self.getPick(traceID))
xaxis = np.arange(min(x)+1, max(x), step)
yaxis = np.arange(min(y)+1, max(y), step)
xaxis = np.arange(min(x), max(x), step)
yaxis = np.arange(min(y), max(y), step)
xgrid, ygrid = np.meshgrid(xaxis, yaxis)
zgrid = griddata((x, y), z, (xgrid, ygrid), method = method)
@ -672,7 +738,7 @@ class SeismicShot(object):
plotmethod[method](*args)
def matshow(self, step = 0.5, method = 'linear', ax = None, plotRec = False, annotations = False):
def matshow(self, ax = None, step = 0.5, method = 'linear', plotRec = True, annotations = True, colorbar = True):
'''
Plots a 2D matrix of the interpolated traveltimes. This needs less performance than plot3dttc
@ -682,27 +748,32 @@ class SeismicShot(object):
:param: method (optional), interpolation method; can be 'linear' (default) or 'cubic'
:type: 'string'
:param: plotRec (optional), plot the receiver positions
:param: plotRec (optional), plot the receiver positions (colored scatter plot, should not be
deactivated because there might be receivers that are not inside the interpolated area)
:type: 'logical'
:param: annotations (optional), displays traceIDs as annotations
:type: 'logical'
'''
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
# plt.interactive('True')
x = []
y = []
z = []
x = []; xcut = []
y = []; ycut = []
z = []; zcut = []
tmin, tmax = self.getCut()
for traceID in self.pick.keys():
if self.getPick(traceID) != None:
if self.getFlag(traceID) != 0:
x.append(self.getRecLoc(traceID)[0])
y.append(self.getRecLoc(traceID)[1])
z.append(self.getPick(traceID))
if self.getFlag(traceID) == 0 and self.getPickIncludeRemoved(traceID) is not None:
xcut.append(self.getRecLoc(traceID)[0])
ycut.append(self.getRecLoc(traceID)[1])
zcut.append(self.getPickIncludeRemoved(traceID))
xaxis = np.arange(min(x)+1, max(x), step)
yaxis = np.arange(min(y)+1, max(y), step)
xaxis = np.arange(min(x), max(x), step)
yaxis = np.arange(min(y), max(y), step)
xgrid, ygrid = np.meshgrid(xaxis, yaxis)
zgrid = griddata((x, y), z, (xgrid, ygrid), method='linear')
@ -710,14 +781,28 @@ class SeismicShot(object):
fig = plt.figure()
ax = plt.axes()
ax.imshow(zgrid, interpolation = 'none', extent = [min(x), max(x), min(y), max(y)])
if annotations == True:
for i, traceID in enumerate(self.pick.keys()):
if shot.picks[traceID] != None:
ax.annotate('%s' % traceID, xy=(x[i], y[i]), fontsize = 'x-small')
ax.matshow(zgrid, extent = [min(x), max(x), min(y), max(y)], origin = 'lower')
plt.text(0.45, 0.9, 'shot: %s' %self.getShotnumber(), transform = ax.transAxes)
sc = ax.scatter(x, y, c = z, s = 30, label = 'picked shots', vmin = tmin, vmax = tmax, linewidths = 1.5)
sccut = ax.scatter(xcut, ycut, c = zcut, s = 30, edgecolor = 'm', label = 'cut out shots', vmin = tmin, vmax = tmax, linewidths = 1.5)
if colorbar == True:
plt.colorbar(sc)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.plot(self.getSrcLoc()[0], self.getSrcLoc()[1],'*k', markersize = 15) # plot source location
if plotRec == True:
ax.plot(x, y, 'k.')
ax.scatter(x, y, c = z, s = 30)
if annotations == True:
for traceID in self.getTraceIDlist():
if self.getFlag(traceID) is not 0:
ax.annotate(' %s' %traceID , xy = (self.getRecLoc(traceID)[0], self.getRecLoc(traceID)[1]),
fontsize = 'x-small', color = 'k')
else:
ax.annotate(' %s' %traceID , xy = (self.getRecLoc(traceID)[0], self.getRecLoc(traceID)[1]),
fontsize = 'x-small', color = 'r')
plt.show()

View File

@ -1,49 +1,191 @@
# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import math
import numpy as np
plt.interactive(True)
class regions(object):
def __init__(self, ax, shot_dict):
'''
A class used for manual inspection and processing of all picks for the user.
Examples:
regions.chooseRectangles():
- lets the user choose several rectangular regions in the plot
regions.plotTracesInRegions():
- creates plots (shot.plot_traces) for all traces in the active regions (i.e. chosen by e.g. chooseRectangles)
regions.setActiveRegionsForDeletion():
- highlights all shots in a the active regions for deletion
regions.deleteMarkedPicks():
- deletes the picks (pick flag set to 0) for all shots set for deletion
regions.deselectSelection(number):
- deselects the region of number = number
'''
def __init__(self, ax, survey):
self.ax = ax
self.shot_dict = shot_dict
self.survey = survey
self.shot_dict = self.survey.getShotDict()
self._x0 = []
self._y0 = []
self._x1 = []
self._y1 = []
self._polyx = []
self._polyy = []
self._allpicks = None
self.shots_found = {}
self.shots_for_deletion = {}
self._generateList()
def _onselect(self, eclick, erelease):
'eclick and erelease are matplotlib events at press and release' #print ' startposition : (%f, %f)' % (eclick.xdata, eclick.ydata)
def _onselect_clicks(self, eclick, erelease):
'eclick and erelease are matplotlib events at press and release'
#print ' startposition : (%f, %f)' % (eclick.xdata, eclick.ydata)
#print ' endposition : (%f, %f)' % (erelease.xdata, erelease.ydata)
print 'region selected x0, y0 = (%3s, %3s), x1, y1 = (%3s, %3s)'%(eclick.xdata, eclick.ydata, erelease.xdata, erelease.ydata)
x0 = min(eclick.xdata, erelease.xdata)
x1 = max(eclick.xdata, erelease.xdata)
y0 = min(eclick.ydata, erelease.ydata)
y1 = max(eclick.ydata, erelease.ydata)
self._x0.append(x0)
self._x1.append(x1)
self._y0.append(y0)
self._y1.append(y1)
self.markCurrentRegion(x0, x1, y0, y1)
shots, numtraces = self.findTracesInShotDict((x0, x1), (y0, y1))
print('Found %d traces in rectangle: %s' %(numtraces, shots))
key = self.getKey()
self.shots_found[key] = {'shots': shots,
'selection': 'rect',
'xvalues': (x0, x1),
'yvalues': (y0, y1)}
self.markRectangle((x0, x1), (y0, y1), key)
def _onselect_verts(self, verts):
x = verts[0][0]
y = verts[0][1]
self._polyx.append(x)
self._polyy.append(y)
self.drawPolyLine()
def _onpress(self, event):
if event.button == 3:
self.disconnectPoly()
def getKey(self):
if self.shots_found.keys() == []:
key = 1
else:
key = max(self.shots_found.keys()) + 1
return key
def drawPolyLine(self):
x = self._polyx
y = self._polyy
if len(x) >= 2 and len(y) >= 2:
plt.plot(x[-2:], y[-2:], 'k')
def drawLastPolyLine(self):
x = self._polyx
y = self._polyy
if len(x) >= 2 and len(y) >= 2:
plt.plot((x[-1], x[0]), (y[-1], y[0]), 'k')
def finishPolygon(self):
self.drawLastPolyLine()
x = self._polyx
y = self._polyy
self._polyx = []; self._polyy = []
key = self.getKey()
self.markPolygon(x, y, key = key)
shots, numtraces = self.findTracesInPoly(x, y)
self.shots_found[key] = {'shots': shots,
'selection': 'poly',
'xvalues': x,
'yvalues': y}
print('Found %d traces in polygon: %s' %(numtraces, shots))
def markPolygon(self, x, y, key = None, color = 'grey', alpha = 0.1, linewidth = 1):
from matplotlib.patches import Polygon
poly = Polygon(np.array(zip(x, y)), color = color, alpha = alpha, lw = linewidth)
self.ax.add_patch(poly)
if key is not None:
self.ax.text((min(x) + (max(x) - min(x)) / 2), (min(y) + (max(y) - min(y)) / 2), str(key))
self.drawFigure()
def disconnectPoly(self):
self.ax.figure.canvas.mpl_disconnect(self._cid)
del self._cid
self.finishPolygon()
self._lasso.disconnect_events()
print('disconnected poly selection\n')
def disconnectRect(self):
self.ax.figure.canvas.mpl_disconnect(self._cid)
del self._cid
self._rectangle.disconnect_events()
print('disconnected rectangle selection\n')
def chooseRectangles(self):
'''
Activates matplotlib widget RectangleSelector.
'''
from matplotlib.widgets import RectangleSelector
print 'Select rectangle is active'
return RectangleSelector(self.ax, self._onselect)
print('Select rectangle is active')
self._cid = self.ax.figure.canvas.mpl_connect('button_press_event', self._onpress)
self._rectangle = RectangleSelector(self.ax, self._onselect_clicks)
return self._rectangle
def _getx0(self):
return self._x0
def choosePolygon(self):
'''
Activates matplotlib widget LassoSelector.
'''
from matplotlib.widgets import LassoSelector
def _getx1(self):
return self._x1
print('Select polygon is active')
self._cid = self.ax.figure.canvas.mpl_connect('button_press_event', self._onpress)
self._lasso = LassoSelector(self.ax, self._onselect_verts)
return self._lasso
def _gety0(self):
return self._y0
def deselectLastSelection(self):
if self.shots_found.keys() == []:
print('No selection found.')
return
key = max(self.shots_found.keys())
self.deselectSelection(key)
def _gety1(self):
return self._y1
def deselectSelection(self, key, color = 'green', alpha = 0.1):
if not key in self.shots_found.keys():
print('No selection found.')
return
if color is not None:
if self.shots_found[key]['selection'] == 'rect':
self.markRectangle(self.shots_found[key]['xvalues'],
self.shots_found[key]['yvalues'],
key = key, color = color, alpha = alpha,
linewidth = 1)
elif self.shots_found[key]['selection'] == 'poly':
self.markPolygon(self.shots_found[key]['xvalues'],
self.shots_found[key]['yvalues'],
key = key, color = color, alpha = alpha,
linewidth = 1)
value = self.shots_found.pop(key)
print('Deselected selection number %d'% key)
return
def _generateList(self):
allpicks = []
for shot in self.shot_dict.values():
for traceID in shot.getTraceIDlist():
allpicks.append((shot.getDistance(traceID), shot.getPickIncludeRemoved(traceID),
shot.getShotnumber(), traceID, shot.getFlag(traceID)))
allpicks.sort()
self._allpicks = allpicks
def getShotDict(self):
return self.shot_dict
@ -51,119 +193,216 @@ class regions(object):
def getShotsForDeletion(self):
return self.shots_for_deletion
def findTracesInShotDict(self, picks = 'normal'):
'''
Returns traces corresponding to a certain area in a plot with all picks over the distances.
'''
print "findTracesInShotDict: Searching for marked traces in the shot dictionary... "
def findTracesInPoly(self, x, y, picks = 'normal', highlight = True):
def dotproduct(v1, v2):
return sum((a*b) for a, b in zip(v1, v2))
for shot in self.shot_dict.values():
whichpicks = {'normal': shot.getPick,
'includeCutOut': shot.getPick_backup}
for index in range(len(self._getx1())):
distancebin = (self._getx0()[index], self._getx1()[index])
pickbin = (self._gety0()[index], self._gety1()[index])
if shot.getTraceIDs4Dist(distancebin = distancebin) is not None:
for traceID in shot.getTraceIDs4Dist(distancebin = distancebin):
if pickbin[0] < whichpicks[picks](traceID) < pickbin[1]:
def getlength(v):
return math.sqrt(dotproduct(v, v))
def getangle(v1, v2):
return np.rad2deg(math.acos(dotproduct(v1, v2) / (getlength(v1) * getlength(v2))))
def insidePoly(x, y, pickX, pickY):
angle = 0
epsilon = 10e-8
for index in range(len(x)):
xval1 = x[index - 1]; yval1 = y[index - 1]
xval2 = x[index]; yval2 = y[index]
angle += getangle([xval1 - pickX, yval1 - pickY], [xval2 - pickX, yval2 - pickY])
if 360 - epsilon <= angle <= 360 + epsilon: ### IMPROVE THAT??
return True
if len(x) == 0 or len(y) == 0:
print('No polygon defined.')
return
shots_found = {}; numtraces = 0
x0 = min(x); x1 = max(x)
y0 = min(y); y1 = max(y)
shots, numtracesrect = self.findTracesInShotDict((x0, x1), (y0, y1), highlight = False)
for shotnumber in shots.keys():
shot = self.shot_dict[shotnumber]
for traceID in shots[shotnumber]:
if shot.getFlag(traceID) is not 0:
pickX = shot.getDistance(traceID)
pickY = shot.getPick(traceID)
if insidePoly(x, y, pickX, pickY):
if not shotnumber in shots_found.keys():
shots_found[shotnumber] = []
shots_found[shotnumber].append(traceID)
if highlight == True:
self.highlightPick(shot, traceID)
if shot.getShotnumber() not in self.shots_found.keys():
self.shots_found[shot.getShotnumber()] = []
if traceID not in self.shots_found[shot.getShotnumber()]:
self.shots_found[shot.getShotnumber()].append(traceID)
self.refreshFigure()
print self.shots_found
numtraces += 1
self.drawFigure()
return shots_found, numtraces
def findTracesInShotDict(self, (x0, x1), (y0, y1), picks = 'normal', highlight = True):
'''
Returns traces corresponding to a certain area in the plot with all picks over the distances.
'''
shots_found = {}; numtraces = 0
if picks == 'normal': pickflag = 0
elif picks == 'includeCutOut': pickflag = None
for line in self._allpicks:
dist, pick, shotnumber, traceID, flag = line
if flag == pickflag: continue ### IMPROVE THAT
if (x0 <= dist <= x1 and y0 <= pick <= y1):
if not shotnumber in shots_found.keys():
shots_found[shotnumber] = []
shots_found[shotnumber].append(traceID)
if highlight == True:
self.highlightPick(self.shot_dict[shotnumber], traceID)
numtraces += 1
self.drawFigure()
return shots_found, numtraces
def highlightPick(self, shot, traceID, annotations = True):
'''
Highlights a single pick for a shot(object)/shotnumber and traceID.
If annotations == True: Displays shotnumber and traceID in the plot.
'''
if type(shot) == int:
shot = self.survey.getShotDict()[shot]
self.ax.scatter(shot.getDistance(traceID), shot.getPick(traceID), s = 50, marker = 'o', facecolors = 'none', edgecolors = 'm', alpha = 1)
if annotations == True:
self.ax.annotate(s = 's%s|t%s'%(shot.getShotnumber(), traceID), xy = (shot.getDistance(traceID), shot.getPick(traceID)), fontsize = 'xx-small')
self.ax.set_ylim(shot.getCut())
def plotTracesInRegion(self):
def highlightAllRegions(self):
'''
Highlights all picks in all active regions.
'''
for key in self.shots_found.keys():
for shotnumber in self.shots_found[key]['shots'].keys():
for traceID in self.shots_found[key]['shots'][shotnumber]:
self.highlightPick(self.shot_dict[shotnumber], traceID)
self.drawFigure()
def plotTracesInRegions(self, keys = 'all', maxfigures = 20):
'''
Plots all traces in the active region or for all specified keys.
:param: keys
:type: int or list
:param: maxfigures, maximum value of figures opened
:type: int
'''
count = 0
maxfigures = 20
# if len(self.shots_found) == 0:
self.findTracesInShotDict()
if keys == 'all':
keys = self.shots_found.keys()
elif type(keys) == int:
keys = [keys]
if len(self.shots_found) > 0:
for shot in self.shot_dict.values():
for shotnumber in self.shots_found:
for key in keys:
for shotnumber in self.shots_found[key]['shots']:
if shot.getShotnumber() == shotnumber:
for traceID in self.shots_found[shotnumber]:
for traceID in self.shots_found[key]['shots'][shotnumber]:
count += 1
if count > maxfigures:
print 'Maximum number of figures (%s) reached. %sth figure was not opened.' %(maxfigures, count)
break
shot.plot_traces(traceID)
else:
print 'No picks yet defined in the regions x = (%s, %s), y = (%s, %s)' %(self._x0, self._x1, self._y0, self._y1)
print('No picks defined in that region(s)')
def plotTracesInRegion_withCutOutTraces(self):
count = 0
maxfigures = 20
# if len(self.shots_found) == 0:
self.findTracesInShotDict(picks = 'includeCutOut')
def setActiveRegionsForDeletion(self):
keys = []
for key in self.shots_found.keys():
keys.append(key)
self.setRegionForDeletion(keys)
if len(self.shots_found) > 0:
for shot in self.shot_dict.values():
for shotnumber in self.shots_found:
if shot.getShotnumber() == shotnumber:
for traceID in self.shots_found[shotnumber]:
count += 1
if count > maxfigures:
print 'Maximum number of figures (%s) reached. %sth figure was not opened.' %(maxfigures, count)
break
shot.plot_traces(traceID)
else:
print 'No picks yet defined in the regions x = (%s, %s), y = (%s, %s)' %(self._x0, self._x1, self._y0, self._y1)
def setRegionForDeletion(self, keys):
if type(keys) == int:
keys = [keys]
def setCurrentRegionsForDeletion(self):
# if len(self.shots_found) == 0:
self.findTracesInShotDict()
for shotnumber in self.shots_found:
for key in keys:
for shotnumber in self.shots_found[key]['shots'].keys():
if not shotnumber in self.shots_for_deletion:
self.shots_for_deletion[shotnumber] = []
for traceID in self.shots_found[shotnumber]:
for traceID in self.shots_found[key]['shots'][shotnumber]:
if not traceID in self.shots_for_deletion[shotnumber]:
self.shots_for_deletion[shotnumber].append(traceID)
self.markAllRegions(color = 'red')
print 'Marked regions for deletion'
self.deselectSelection(key, color = 'red', alpha = 0.2)
def markAllRegions(self, color = 'grey'):
print 'Set region(s) %s for deletion'%keys
def markAllActiveRegions(self):
for key in self.shots_found.keys():
if self.shots_found[key]['selection'] == 'rect':
self.markRectangle(self.shots_found[key]['xvalues'],
self.shots_found[key]['yvalues'], key = key)
if self.shots_found[key]['selection'] == 'poly':
self.markPolygon(self.shots_found[key]['xvalues'],
self.shots_found[key]['yvalues'], key = key)
def markRectangle(self, (x0, x1), (y0, y1), key = None, color = 'grey', alpha = 0.1, linewidth = 1):
'''
Mark a rectangular region on the axes.
'''
from matplotlib.patches import Rectangle
for index in range(len(self._getx0())):
x0 = self._getx0()[index]
y0 = self._gety0()[index]
x1 = self._getx1()[index]
y1 = self._gety1()[index]
self.ax.add_patch(Rectangle((x0, y0), (x1 - x0), (y1 - y0),
alpha = alpha, facecolor = color, linewidth = linewidth))
if key is not None:
self.ax.text((x0 + (x1 - x0) / 2), (y0 + (y1 - y0) / 2), str(key))
self.drawFigure()
self.ax.add_patch(Rectangle((x0, y0), (x1 - x0), (y1 - y0), alpha=0.5, facecolor = color))
self.refreshFigure()
def refreshFigure(self):
print('Refreshing figure...')
self.ax.clear()
self.ax = self.survey.plotAllPicks(ax = self.ax, refreshPlot = True)
self.markAllActiveRegions()
self.drawFigure()
print('Done!')
def markCurrentRegion(self, x0, x1, y0, y1, color = 'grey'):
from matplotlib.patches import Rectangle
def clearShotsForDeletion(self):
'''
Clears the list of shots marked for deletion.
'''
self.shots_for_deletion = {}
print('Cleared all shots that were set for deletion.')
self.ax.add_patch(Rectangle((x0, y0), (x1 - x0), (y1 - y0), alpha=0.1, facecolor = color))
self.refreshFigure()
def getShotsForDeletion(self):
return self.shots_for_deletion
def deleteMarkedPicks(self):
'''
Deletes all shots set for deletion.
'''
if len(self.getShotsForDeletion()) is 0:
print('No shots set for deletion.')
return
for shot in self.getShotDict().values():
for shotnumber in self.getShotsForDeletion():
if shot.getShotnumber() == shotnumber:
for traceID in self.getShotsForDeletion()[shotnumber]:
shot.removePick(traceID)
print "Deleted the pick for traceID %s on shot number %s" %(traceID, shotnumber)
self.shots_for_deletion = {} # clear dictionary
def highlightPicksForShot(self, shot, annotations = False):
for traceID in shot.getTraceIDlist():
if shot.getPick(traceID) is not None:
self.highlightPick(shot, traceID, annotations)
self.clearShotsForDeletion()
self.refreshFigure()
def refreshFigure(self):
def highlightPicksForShot(self, shot, annotations = False):
'''
Highlight all picks for a given shot.
'''
if type(shot) is int:
shot = self.survey.getShotDict()[shotnumber]
for traceID in shot.getTraceIDlist():
if shot.getFlag(traceID) is not 0:
self.highlightPick(shot, traceID, annotations)
self.drawFigure()
def drawFigure(self):
plt.draw()

View File

@ -73,7 +73,7 @@ def fitSNR4dist(shot_dict, shiftdist = 5):
for traceID in shot.getTraceIDlist():
if shot.getSNR(traceID)[0] >= 1:
dists.append(shot.getDistance(traceID))
picks.append(shot.getPick_backup(traceID))
picks.append(shot.getPickIncludeRemoved(traceID))
snrs.append(shot.getSNR(traceID)[0])
snr_sqrt_inv.append(1/np.sqrt(shot.getSNR(traceID)[0]))
fit = np.polyfit(dists, snr_sqrt_inv, 1)

View File

@ -25,7 +25,7 @@ class CharacteristicFunction(object):
'''
SuperClass for different types of characteristic functions.
'''
def __init__(self, data, cut, t2=None, order=None, t1=None, fnoise=None):
def __init__(self, data, cut, t2=None, order=None, t1=None, fnoise=None, stealthMode=False):
'''
Initialize data type object with information from the original
Seismogram.
@ -62,6 +62,7 @@ class CharacteristicFunction(object):
self.calcCF(self.getDataArray())
self.arpara = np.array([])
self.xpred = np.array([])
self._stealthMode = stealthMode
def __str__(self):
return '''\n\t{name} object:\n
@ -135,6 +136,9 @@ class CharacteristicFunction(object):
def getXCF(self):
return self.xcf
def _getStealthMode(self):
return self._stealthMode()
def getDataArray(self, cut=None):
'''
If cut times are given, time series is cut from cut[0] (start time)
@ -219,7 +223,8 @@ class AICcf(CharacteristicFunction):
def calcCF(self, data):
#print 'Calculating AIC ...' ## MP MP output suppressed
#if self._getStealthMode() is False:
# print 'Calculating AIC ...'
x = self.getDataArray()
xnp = x[0].data
nn = np.isnan(xnp)
@ -257,11 +262,13 @@ class HOScf(CharacteristicFunction):
if len(nn) > 1:
xnp[nn] = 0
if self.getOrder() == 3: # this is skewness
print 'Calculating skewness ...'
#if self._getStealthMode() is False:
# print 'Calculating skewness ...'
y = np.power(xnp, 3)
y1 = np.power(xnp, 2)
elif self.getOrder() == 4: # this is kurtosis
#print 'Calculating kurtosis ...' ## MP MP output suppressed
#if self._getStealthMode() is False:
# print 'Calculating kurtosis ...'
y = np.power(xnp, 4)
y1 = np.power(xnp, 2)

View File

@ -15,7 +15,7 @@ from obspy.core import Stream, UTCDateTime
import warnings
def earllatepicker(X, nfac, TSNR, Pick1, iplot=None):
def earllatepicker(X, nfac, TSNR, Pick1, iplot=None, stealthMode = False):
'''
Function to derive earliest and latest possible pick after Diehl & Kissling (2009)
as reasonable uncertainties. Latest possible pick is based on noise level,
@ -44,7 +44,8 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None):
LPick = None
EPick = None
PickError = None
#print 'earllatepicker: Get earliest and latest possible pick relative to most likely pick ...'
if stealthMode is False:
print 'earllatepicker: Get earliest and latest possible pick relative to most likely pick ...'
x = X[0].data
t = np.arange(0, X[0].stats.npts / X[0].stats.sampling_rate,
@ -75,6 +76,7 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None):
# if EPick stays NaN the signal window size will be doubled
while np.isnan(EPick):
if count > 0:
if stealthMode is False:
print("\nearllatepicker: Doubled signal window size %s time(s) "
"because of NaN for earliest pick." %count)
isigDoubleWinStart = pis[-1] + 1