diff --git a/QtPyLoT.py b/QtPyLoT.py index cae007c4..2ba21c67 100755 --- a/QtPyLoT.py +++ b/QtPyLoT.py @@ -273,8 +273,7 @@ class MainWindow(QMainWindow): slot=self.autoPick, shortcut='Alt+Ctrl+A', icon=auto_icon, tip='Automatically pick' ' the entire dataset' - ' displayed!', - checkable=False) + ' displayed!') autoPickToolBar = self.addToolBar("autoPyLoT") autoPickActions = (auto_pick,) diff --git a/autoPyLoT.py b/autoPyLoT.py index 767881af..af61bcb1 100755 --- a/autoPyLoT.py +++ b/autoPyLoT.py @@ -8,12 +8,11 @@ import glob import matplotlib.pyplot as plt from obspy.core import read -from pylot.core.util import _getVersionString from pylot.core.read.data import Data from pylot.core.read.inputs import AutoPickParameter from pylot.core.util.structure import DATASTRUCTURE from pylot.core.pick.autopick import autopickevent -from pylot.core.pick.utils import writephases +from pylot.core.util.version import get_git_version as _getVersionString __version__ = _getVersionString() @@ -86,9 +85,6 @@ def autoPyLoT(inputfile): # !automated picking starts here! picks = autopickevent(wfdat, parameter) - # write phases to NLLoc-phase file - writephases(wd_checked_onsets, 'NLLoc', phasefile) - print '------------------------------------------' print '-----Finished event %s!-----' % event print '------------------------------------------' @@ -104,9 +100,6 @@ def autoPyLoT(inputfile): # !automated picking starts here! picks = autopickevent(wfdat, parameter) - # write phases to NLLoc-phase file - writephases(wd_checked_onsets, 'NLLoc', phasefile) - print '------------------------------------------' print '-------Finished event %s!-------' % parameter.getParam('eventID') print '------------------------------------------' diff --git a/icons.qrc b/icons.qrc index efd43234..010efb1b 100644 --- a/icons.qrc +++ b/icons.qrc @@ -1,10 +1,10 @@ - icons/pylot.ico - icons/pylot.png + icons/pylot.ico + icons/pylot.png icons/printer.png icons/delete.png - icons/key_E.png + icons/key_E.png icons/key_N.png icons/key_P.png icons/key_Q.png @@ -14,7 +14,7 @@ icons/key_U.png icons/key_V.png icons/key_W.png - icons/key_Z.png + icons/key_Z.png icons/filter.png icons/sync.png icons/zoom_0.png diff --git a/pylot/core/__init__.py b/pylot/core/__init__.py index e69de29b..40a96afc 100755 --- a/pylot/core/__init__.py +++ b/pylot/core/__init__.py @@ -0,0 +1 @@ +# -*- coding: utf-8 -*- diff --git a/pylot/core/active/__init__.py b/pylot/core/active/__init__.py index ccbcc844..9ef8ae9f 100644 --- a/pylot/core/active/__init__.py +++ b/pylot/core/active/__init__.py @@ -1 +1,2 @@ +# -*- coding: utf-8 -*- __author__ = 'sebastianw' diff --git a/pylot/core/active/activeSeismoPick.py b/pylot/core/active/activeSeismoPick.py index 50c1ddad..3dc97353 100644 --- a/pylot/core/active/activeSeismoPick.py +++ b/pylot/core/active/activeSeismoPick.py @@ -1,3 +1,4 @@ +# -*- coding: utf-8 -*- import sys import numpy as np from pylot.core.active import seismicshot @@ -14,11 +15,12 @@ class Survey(object): self._sourcefile = sourcefile self._obsdir = path self._generateSurvey() - if useDefaultParas == True: + if useDefaultParas == True: self.setParametersForShots() self._removeAllEmptyTraces() self._updateShots() - + self.setArtificialPick(0, 0) + def _generateSurvey(self): from obspy.core import read @@ -27,23 +29,28 @@ class Survey(object): for shotnumber in shotlist: # loop over data files # generate filenames and read manual picks to a list obsfile = self._obsdir + str(shotnumber) + '_pickle.dat' - - if not obsfile in shot_dict.keys(): + if obsfile not in shot_dict.keys(): shot_dict[shotnumber] = [] 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) @@ -65,7 +72,7 @@ class Survey(object): if removed is not None: if count == 0: outfile = open(filename, 'w') count += 1 - outfile.writelines('shot: %s, removed empty traces: %s\n' + outfile.writelines('shot: %s, removed empty traces: %s\n' %(shot.getShotnumber(), removed)) print ("\nremoveEmptyTraces: Finished! Removed %d traces" %count) if count > 0: @@ -83,7 +90,7 @@ class Survey(object): count += 1 countTraces += len(del_traceIDs) outfile.writelines("shot: %s, removed traceID(s) %s because " - "they were not found in the corresponding stream\n" + "they were not found in the corresponding stream\n" %(shot.getShotnumber(), del_traceIDs)) print ("\nupdateShots: Finished! Updated %d shots and removed " @@ -120,14 +127,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 @@ -137,6 +143,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) @@ -144,8 +166,8 @@ class Survey(object): def countAllTraces(self): numtraces = 0 - for line in self.getShotlist(): - for line in self.getReceiverlist(): + for shot in self.getShotlist(): + for rec in self.getReceiverlist(): ### shot.getReceiverlist etc. numtraces += 1 return numtraces @@ -180,7 +202,7 @@ class Survey(object): def getReceiverfile(self): return self._recfile - + def getPath(self): return self._obsdir @@ -194,7 +216,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, @@ -228,14 +250,14 @@ class Survey(object): (x, y, z) = shot.getSrcLoc() # getSrcLoc returns (x, y, z) srcfile.writelines('%10s %10s %10s\n' %(getAngle(y), getAngle(x), (-1)*z)) # lat, lon, depth LatAll.append(getAngle(y)); LonAll.append(getAngle(x)); DepthAll.append((-1)*z) - srcfile.writelines('%10s\n' %1) # + srcfile.writelines('%10s\n' %1) # srcfile.writelines('%10s %10s %10s\n' %(1, 1, ttfilename)) ttfile = open(directory + '/' + ttfilename, 'w') traceIDlist = shot.getTraceIDlist() 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) @@ -243,7 +265,7 @@ class Survey(object): LatAll.append(getAngle(y)); LonAll.append(getAngle(x)); DepthAll.append((-1)*z) count += 1 ttfile.close() - srcfile.close() + srcfile.close() print 'Wrote output for %s traces' %count print 'WARNING: output generated for FMTOMO-obsdata. Obsdata seems to take Lat, Lon, Depth and creates output for FMTOMO as Depth, Lat, Lon' print 'Dimensions of the seismic Array, transformed for FMTOMO, are Depth(%s, %s), Lat(%s, %s), Lon(%s, %s)'%( @@ -252,14 +274,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, cbar = 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) @@ -267,44 +362,63 @@ 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) - ax.legend() + color = {'log10SNR': snrlog, + 'pickerror': pickerror, + 'spe': spe} + self.color = color + if refreshPlot is False: + ax, cbar = self.createPlot(dist, pick, color[colorByVal], label='%s' % colorByVal) + region = regions(ax, cbar, self) + ax.legend() + return (ax, region) + if refreshPlot is True: + ax, cbar = self.createPlot(dist, pick, color[colorByVal], label='%s' % colorByVal, ax=ax, cbar=cbar) + ax.legend() + return ax - return ax, region - - def createPlot(self, dist, pick, inkByVal, label): + def createPlot(self, dist, pick, inkByVal, label, ax = None, cbar = None): import matplotlib.pyplot as plt plt.interactive(True) cm = plt.cm.jet - - fig = plt.figure() - ax = fig.add_subplot(111) - fig = ax.scatter(dist, pick, cmap = cm, c = inkByVal, s = 5, edgecolors = 'none', label = label) - cbar = plt.colorbar(fig, fraction = 0.05) - cbar.set_label(label) - plt.title('Plot of all Picks') - plt.xlabel('Distance [m]') - plt.ylabel('Time [s]') - - return ax + if ax is None: + print('Generating new plot...') + fig = plt.figure() + ax = fig.add_subplot(111) + sc = ax.scatter(dist, pick, cmap=cm, c=inkByVal, s=5, edgecolors='none', label=label) + cbar = plt.colorbar(sc, fraction=0.05) + cbar.set_label(label) + ax.set_xlabel('Distance [m]') + ax.set_ylabel('Time [s]') + ax.text(0.5, 0.95, 'Plot of all picks', transform=ax.transAxes, horizontalalignment='center') + else: + sc = ax.scatter(dist, pick, cmap=cm, c=inkByVal, s=5, edgecolors='none', label=label) + cbar = plt.colorbar(sc, cax=cbar.ax) + cbar.set_label(label) + ax.set_xlabel('Distance [m]') + ax.set_ylabel('Time [s]') + ax.text(0.5, 0.95, 'Plot of all picks', transform=ax.transAxes, horizontalalignment='center') + return (ax, cbar) def _update_progress(self, shotname, tend, progress): - sys.stdout.write("Working on shot %s. ETC is %02d:%02d:%02d [%2.2f %%]\r" - %(shotname, tend.hour, tend.minute, tend.second, progress)) + sys.stdout.write('Working on shot %s. ETC is %02d:%02d:%02d [%2.2f %%]\r' % (shotname, + tend.hour, + tend.minute, + tend.second, + progress)) sys.stdout.flush() def saveSurvey(self, filename = 'survey.pickle'): @@ -313,8 +427,8 @@ class Survey(object): cPickle.dump(self, outfile, -1) print('saved Survey to file %s'%(filename)) - - @staticmethod + + @staticmethod def from_pickle(filename): import cPickle infile = open(filename, 'rb') diff --git a/pylot/core/active/fmtomo2vtk.py b/pylot/core/active/fmtomo2vtk.py index a4edf32d..95b90705 100644 --- a/pylot/core/active/fmtomo2vtk.py +++ b/pylot/core/active/fmtomo2vtk.py @@ -1,3 +1,4 @@ +# -*- coding: utf-8 -*- import numpy as np def vgrids2VTK(inputfile = 'vgrids.in', outputfile = 'vgrids.vtk'): @@ -73,7 +74,7 @@ def vgrids2VTK(inputfile = 'vgrids.in', outputfile = 'vgrids.vtk'): dX = getDistance(np.rad2deg(dPhi)) dY = getDistance(np.rad2deg(dTheta)) - + nPoints = nX * nY * nZ dZ = dR @@ -114,7 +115,7 @@ def rays2VTK(fnin, fdirout = './vtk_files/', nthPoint = 50): R = 6371. distance = angle / 180 * (PI * R) return distance - + infile = open(fnin, 'r') R = 6371 rays = {} @@ -124,10 +125,15 @@ def rays2VTK(fnin, fdirout = './vtk_files/', nthPoint = 50): ### NOTE: rays.dat seems to be in km and radians while True: - raynumber += 1 + 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] = {} diff --git a/pylot/core/active/picking_script.py b/pylot/core/active/picking_script.py deleted file mode 100644 index fbb3ef62..00000000 --- a/pylot/core/active/picking_script.py +++ /dev/null @@ -1,106 +0,0 @@ -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 diff --git a/pylot/core/active/seismicArrayPreparation.py b/pylot/core/active/seismicArrayPreparation.py index 65226423..7e7dd09a 100644 --- a/pylot/core/active/seismicArrayPreparation.py +++ b/pylot/core/active/seismicArrayPreparation.py @@ -1,3 +1,4 @@ +# -*- coding: utf-8 -*- import sys import numpy as np from scipy.interpolate import griddata @@ -28,12 +29,12 @@ class SeisArray(object): def _generateReceiverlines(self): ''' - Connects the traceIDs to the lineIDs + 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]) + 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) @@ -43,16 +44,16 @@ class SeisArray(object): 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]) + 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]) + traceID = int(line.split()[0]) gphoneNum = float(line.split()[2]) self._geophoneNumbers[traceID] = gphoneNum @@ -93,7 +94,7 @@ class SeisArray(object): return self._geophoneNumbers[traceID] def getMeasuredReceivers(self): - return self._measuredReceivers + return self._measuredReceivers def getMeasuredTopo(self): return self._measuredTopo @@ -139,11 +140,11 @@ class SeisArray(object): if self._getReceiverValue(traceID1, coordinate) < self._getReceiverValue(traceID2, coordinate): direction = +1 return direction - if self._getReceiverValue(traceID1, coordinate) > self._getReceiverValue(traceID2, coordinate): + 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 @@ -186,7 +187,7 @@ class SeisArray(object): x = float(line[1]) y = float(line[2]) z = float(line[3]) - self._measuredTopo[pointID] = (x, y, z) + self._measuredTopo[pointID] = (x, y, z) def addSourceLocations(self, filename): ''' @@ -202,7 +203,7 @@ class SeisArray(object): x = float(line[1]) y = float(line[2]) z = float(line[3]) - self._sourceLocs[pointID] = (x, y, z) + self._sourceLocs[pointID] = (x, y, z) def interpZcoords4rec(self, method = 'linear'): ''' @@ -239,9 +240,9 @@ class SeisArray(object): ''' x = []; y = []; z = [] for traceID in self.getMeasuredReceivers().keys(): - x.append(self.getMeasuredReceivers()[traceID][0]) + x.append(self.getMeasuredReceivers()[traceID][0]) y.append(self.getMeasuredReceivers()[traceID][1]) - z.append(self.getMeasuredReceivers()[traceID][2]) + z.append(self.getMeasuredReceivers()[traceID][2]) return x, y, z def getMeasuredTopoLists(self): @@ -250,9 +251,9 @@ class SeisArray(object): ''' x = []; y = []; z = [] for pointID in self.getMeasuredTopo().keys(): - x.append(self.getMeasuredTopo()[pointID][0]) + x.append(self.getMeasuredTopo()[pointID][0]) y.append(self.getMeasuredTopo()[pointID][1]) - z.append(self.getMeasuredTopo()[pointID][2]) + z.append(self.getMeasuredTopo()[pointID][2]) return x, y, z def getSourceLocsLists(self): @@ -261,9 +262,9 @@ class SeisArray(object): ''' x = []; y = []; z = [] for pointID in self.getSourceLocations().keys(): - x.append(self.getSourceLocations()[pointID][0]) + x.append(self.getSourceLocations()[pointID][0]) y.append(self.getSourceLocations()[pointID][1]) - z.append(self.getSourceLocations()[pointID][2]) + z.append(self.getSourceLocations()[pointID][2]) return x, y, z def getAllMeasuredPointsLists(self): @@ -289,7 +290,7 @@ class SeisArray(object): 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. @@ -317,7 +318,7 @@ class SeisArray(object): :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 @@ -356,9 +357,9 @@ class SeisArray(object): progress = float(count) / float(nTotal) * 100 self._update_progress(progress) - if filename is not None: + if filename is not None: outfile.writelines('%10s\n'%(z + elevation)) - + return surface def generateVgrid(self, nTheta = 80, nPhi = 80, nR = 120, @@ -415,7 +416,7 @@ class SeisArray(object): 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 @@ -455,7 +456,7 @@ class SeisArray(object): progress = float(count) / float(nTotal) * 100 self._update_progress(progress) - + outfile.close() def exportAll(self, filename = 'interpolated_receivers.out'): @@ -463,7 +464,7 @@ class SeisArray(object): count = 0 for traceID in self.getReceiverCoordinates().keys(): count += 1 - x, y, z = self.getReceiverCoordinates()[traceID] + 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() @@ -472,15 +473,15 @@ class SeisArray(object): import matplotlib.pyplot as plt plt.interactive(True) plt.figure() - xmt, ymt, zmt = self.getMeasuredTopoLists() + xmt, ymt, zmt = self.getMeasuredTopoLists() xsc, ysc, zsc = self.getSourceLocsLists() - xmr, ymr, zmr = self.getMeasuredReceiverLists() - xrc, yrc, zrc = self.getReceiverLists() + 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: + 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') @@ -490,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 @@ -501,9 +506,9 @@ class SeisArray(object): fig = plt.figure() ax = plt.axes(projection = '3d') - xmt, ymt, zmt = self.getMeasuredTopoLists() - xmr, ymr, zmr = self.getMeasuredReceiverLists() - xin, yin, zin = self.getReceiverLists() + 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') @@ -512,8 +517,8 @@ class SeisArray(object): ax.legend() return ax - - + + def plotSurface3D(self, ax = None, step = 0.5, method = 'linear'): from matplotlib import cm import matplotlib.pyplot as plt @@ -657,8 +662,8 @@ class SeisArray(object): cPickle.dump(self, outfile, -1) print('saved SeisArray to file %s'%(filename)) - - @staticmethod + + @staticmethod def from_pickle(filename): import cPickle infile = open(filename, 'rb') diff --git a/pylot/core/active/seismicArrayPreperation.py b/pylot/core/active/seismicArrayPreperation.py deleted file mode 100644 index ae1b42ec..00000000 --- a/pylot/core/active/seismicArrayPreperation.py +++ /dev/null @@ -1,665 +0,0 @@ -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) diff --git a/pylot/core/active/seismicshot.py b/pylot/core/active/seismicshot.py index ecb602c2..07cd2b43 100644 --- a/pylot/core/active/seismicshot.py +++ b/pylot/core/active/seismicshot.py @@ -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 = {} @@ -44,14 +43,14 @@ class SeismicShot(object): removed = [] for i in range(0, len(coordlist)): traceIDs.append(int(coordlist[i].split()[0])) - + for trace in self.traces: try: traceIDs.index(int(trace.stats.channel)) except: - self.traces.remove(trace) + self.traces.remove(trace) removed.append(int(trace.stats.channel)) - + if len(removed) > 0: return removed @@ -59,7 +58,7 @@ class SeismicShot(object): for trace in self.traces: if traceID == trace.stats.channel: 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: @@ -133,24 +132,34 @@ 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)) - if np.isnan(pickerror) == True: + if np.isnan(pickerror) == True: print("SPE is NaN for shot %s, traceID %s"%(self.getShotnumber(), traceID)) - return pickerror - + return pickerror + def getStreamTraceIDs(self): traceIDs = [] for trace in self.traces: @@ -172,15 +181,15 @@ class SeismicShot(object): def getPickwindow(self, traceID): try: - self.pickwindow[traceID] + self.pickwindow[traceID] except KeyError as e: print('no pickwindow for trace %s, set to %s' % (traceID, self.getCut())) self.setPickwindow(traceID, self.getCut()) return self.pickwindow[traceID] - + def getSNR(self, traceID): return self.snr[traceID] - + def getSNRthreshold(self, traceID): return self.snrthreshold[traceID] @@ -218,7 +227,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): ''' @@ -241,7 +250,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 ? ########## ''' @@ -254,19 +263,19 @@ class SeismicShot(object): traces = [trace for trace in self.traces if int(trace.stats.channel) == traceID] if len(traces) == 1: return Stream(traces) - else: - self.setPick(traceID, None) - print('Warning: ambigious or empty traceID: %s' % traceID) - + self.setPick(traceID, None) + print 'Warning: ambigious or empty traceID: %s' % traceID + #raise ValueError('ambigious or empty traceID: %s' % traceID) - + def pickTraces(self, traceID, windowsize, folm = 0.6, HosAic = 'hos'): ########## input variables ########## + # LOCALMAX NOT IMPLEMENTED! ''' Intitiate picking for a trace. :param: traceID :type: int - + :param: cutwindow (equals HOScf 'cut' variable) :type: tuple @@ -289,29 +298,28 @@ class SeismicShot(object): aiccf = self.getAICcf(traceID) self.timeArray[traceID] = hoscf.getTimeArray() - aiccftime, hoscftime = self.threshold(hoscf, aiccf, windowsize, self.getPickwindow(traceID), folm) - setHosAic = {'hos': hoscftime, '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), - nfac, (tnoise, tgap, tsignal), - self.getPick(traceID)) + (self.pick[traceID]['epp'], self.pick[traceID]['lpp'], + self.pick[traceID]['spe']) = earllatepicker(self.getSingleStream(traceID), + nfac, (tnoise, tgap, tsignal), + self.getPickIncludeRemoved(traceID), + stealthMode = True) def threshold(self, hoscf, aiccf, windowsize, pickwindow, folm = 0.6): ''' - Threshold picker, using the local maximum in a pickwindow to find the time at + Threshold picker, using the local maximum in a pickwindow to find the time at which a fraction of the local maximum is reached for the first time. - + :param: hoscf, Higher Order Statistics Characteristic Function :type: 'Characteristic Function' @@ -337,34 +345,34 @@ class SeismicShot(object): threshold = folm * max(hoscflist[leftb : rightb]) # combination of local maximum and threshold m = leftb - + while hoscflist[m] < threshold: m += 1 - + hoscftime = list(hoscf.getTimeArray())[m] lb = max(0, m - windowsize[0]) # if window exceeds t = 0 aiccfcut = list(aiccf.getCF())[lb : m + windowsize[1]] n = aiccfcut.index(min(aiccfcut)) - + m = lb + n - + aiccftime = list(hoscf.getTimeArray())[m] - + return aiccftime, hoscftime def getDistance(self, traceID): ''' Returns the distance of the receiver with the ID == traceID to the source location (shot location). Uses getSrcLoc and getRecLoc. - + :param: traceID :type: int ''' shotX, shotY, shotZ = self.getSrcLoc() recX, recY, recZ = self.getRecLoc(traceID) dist = np.sqrt((shotX-recX)**2 + (shotY-recY)**2 + (shotZ-recZ)**2) - + if np.isnan(dist) == True: raise ValueError("Distance is NaN for traceID %s" %traceID) @@ -375,7 +383,7 @@ class SeismicShot(object): ''' Returns the location (x, y, z) of the receiver with the ID == traceID. RECEIVEIVER FILE MUST BE SET FIRST, TO BE IMPROVED. - + :param: traceID :type: int ''' @@ -389,7 +397,7 @@ class SeismicShot(object): y = coordlist[i].split()[2] z = coordlist[i].split()[3] return float(x), float(y), float(z) - + #print "WARNING: traceID %s not found" % traceID raise ValueError("traceID %s not found" % traceID) #return float(self.getSingleStream(traceID)[0].stats.seg2['RECEIVER_LOCATION']) @@ -412,7 +420,7 @@ class SeismicShot(object): ''' Returns the traceID(s) for a certain distance between source and receiver. Used for 2D Tomography. TO BE IMPROVED. - + :param: distance :type: real @@ -428,7 +436,7 @@ class SeismicShot(object): if self.getDistance(traceID) == distance: traceID_list.append(traceID) if distancebin[0] >= 0 and distancebin[1] > 0: - if self.getDistance(traceID) > distancebin[0] and self.getDistance(traceID) < distancebin[1]: + if distancebin[0] < self.getDistance(traceID) < distancebin[1]: traceID_list.append(traceID) if len(traceID_list) > 0: @@ -437,7 +445,7 @@ class SeismicShot(object): 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 @@ -452,21 +460,31 @@ class SeismicShot(object): 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) - + # 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 # else: - # raise KeyError('pick to be set more than once for traceID %s' % traceID) + # raise KeyError('pick to be set more than once for traceID %s' % traceID) # def readParameter(self, parfile): # 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 @@ -474,12 +492,13 @@ class SeismicShot(object): def setSNR(self, traceID): ########## FORCED HOS PICK ########## ''' Gets the SNR using pylot and then sets the SNR for the traceID. - + :param: traceID :type: int :param: (tnoise, tgap, tsignal), as used in pylot SNR ''' + from pylot.core.pick.utils import getSNR tgap = self.getTgap() @@ -509,7 +528,7 @@ class SeismicShot(object): # def plot2dttc(self, dist_med = 0): ########## 2D ########## # ''' # Function to plot the traveltime curve for automated picks (AIC & HOS) of a shot. 2d only! - + # :param: dist_med (optional) # :type: 'dictionary' # ''' @@ -517,7 +536,7 @@ class SeismicShot(object): # plt.interactive('True') # aictimearray = [] # hostimearray = [] - # if dist_med is not 0: + # if dist_med is not 0: # dist_medarray = [] # i = 1 @@ -570,42 +589,89 @@ 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 traceID not 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 - + stream = self.getSingleStream(traceID) stime = getGlobalTimes(stream)[0] - timeaxis = prepTimeAxis(stime, stream[0]) + timeaxis = prepTimeAxis(stime, stream[0]) timeaxis -= stime - - plt.interactive('True') + ax = self.traces4plot[traceID]['ax1'] + + 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.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') + ax.legend() + + def _drawCFs(self, traceID, folm, refresh = False): hoscf = self.getHOScf(traceID) aiccf = self.getAICcf(traceID) + ax = self.traces4plot[traceID]['ax2'] - 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)], - [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)], + 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): ''' Plots a 3D 'traveltime cone' as surface plot by interpolating on a regular grid over the traveltimes, not yet regarding the vertical offset of the receivers. @@ -622,7 +688,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 @@ -631,20 +696,20 @@ 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) if ax == None: fig = plt.figure() ax = plt.axes(projection = '3d') - + xsrc, ysrc, zsrc = self.getSrcLoc() if contour == True: @@ -656,13 +721,13 @@ class SeismicShot(object): if plotpicks == True: ax.plot(x, y, z, 'k.') - + def plotttc(self, method, *args): plotmethod = {'2d': self.plot2dttc, '3d': self.plot3dttc} - + 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 @@ -672,27 +737,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') @@ -700,14 +770,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() + + diff --git a/pylot/core/active/surveyPlotTools.py b/pylot/core/active/surveyPlotTools.py index 2f7b2a0a..9a63c21e 100644 --- a/pylot/core/active/surveyPlotTools.py +++ b/pylot/core/active/surveyPlotTools.py @@ -1,48 +1,69 @@ +# -*- 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.plotTracesInActiveRegions(): + - creates plots (shot.plot_traces) for all traces in the active regions (i.e. chosen by e.g. chooseRectangles) + + regions.setAllActiveRegionsForDeletion(): + - highlights all shots in a the active regions for deletion + + regions.deleteAllMarkedPicks(): + - 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, cbar, survey): self.ax = ax - self.shot_dict = shot_dict + self.cbar = cbar + self.cbv = 'log10SNR' + self._xlim0 = self.ax.get_xlim() + self._ylim0 = self.ax.get_ylim() + self._xlim = self.ax.get_xlim() + self._ylim = self.ax.get_ylim() + self.survey = survey + self.shot_dict = self.survey.getShotDict() self._x0 = [] self._y0 = [] self._x1 = [] self._y1 = [] + self._polyx = [] + self._polyy = [] + self.buttons = {} + self._allpicks = None self.shots_found = {} self.shots_for_deletion = {} + self._generateList() + self._addButtons() + self.addTextfield() + self.drawFigure() - def _onselect(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) + 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))) - def chooseRectangles(self): - from matplotlib.widgets import RectangleSelector - - print 'Select rectangle is active' - return RectangleSelector(self.ax, self._onselect) - - def _getx0(self): - return self._x0 - - def _getx1(self): - return self._x1 - - def _gety0(self): - return self._y0 - - def _gety1(self): - return self._y1 + allpicks.sort() + self._allpicks = allpicks def getShotDict(self): return self.shot_dict @@ -50,119 +71,439 @@ 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 _onselect_clicks(self, eclick, erelease): + '''eclick and erelease are matplotlib events at press and release''' + 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) - 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]: + shots, numtraces = self.findTracesInShotDict((x0, x1), (y0, y1)) + self.printOutput('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) + self.disconnectRect() + + 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() + self.printOutput('Disconnected polygon selection') + + def addTextfield(self, xpos = 0, ypos = 0.95, width = 1, height = 0.03): + self.axtext = self.ax.figure.add_axes([xpos, + ypos, + width, + height]) + self.axtext.xaxis.set_visible(False) + self.axtext.yaxis.set_visible(False) + + def writeInTextfield(self, text = None): + self.setXYlim(self.ax.get_xlim(), self.ax.get_ylim()) + self.axtext.clear() + self.axtext.text(0.01, 0.5, text, verticalalignment='center', horizontalalignment='left') + self.drawFigure() + + def _addButtons(self): + xpos1 = 0.13 + xpos2 = 0.6 + dx = 0.06 + self.addButton('Rect', self.chooseRectangles, xpos=xpos1, color='white') + self.addButton('Poly', self.choosePolygon, xpos=xpos1 + dx, color='white') + self.addButton('Plot', self.plotTracesInActiveRegions, xpos=xpos1 + 2 * dx, color='yellow') + self.addButton('SNR', self.refreshLog10SNR, xpos=xpos1 + 3 * dx, color='cyan') + self.addButton('PE', self.refreshPickerror, xpos=xpos1 + 4 * dx, color='cyan') + self.addButton('SPE', self.refreshSPE, xpos=xpos1 + 5 * dx, color='cyan') + self.addButton('DesLst', self.deselectLastSelection, xpos=xpos2 + dx, color='green') + self.addButton('SelAll', self.setAllActiveRegionsForDeletion, xpos=xpos2 + 2 * dx) + self.addButton('DelAll', self.deleteAllMarkedPicks, xpos=xpos2 + 3 * dx, color='red') + + def addButton(self, name, action, xpos, ypos = 0.91, color = None): + from matplotlib.widgets import Button + self.buttons[name] = {'ax': None, + 'button': None, + 'action': action, + 'xpos': xpos} + ax = self.ax.figure.add_axes([xpos, + ypos, + 0.05, + 0.03]) + button = Button(ax, name, color=color, hovercolor='grey') + button.on_clicked(action) + self.buttons[name]['ax'] = ax + self.buttons[name]['button'] = button + self.buttons[name]['xpos'] = xpos + + def getKey(self): + if self.shots_found.keys() == []: + key = 1 + else: + key = max(self.shots_found.keys()) + 1 + return key + + def drawPolyLine(self): + self.setXYlim(self.ax.get_xlim(), self.ax.get_ylim()) + x = self._polyx + y = self._polyy + if len(x) >= 2 and len(y) >= 2: + self.ax.plot(x[-2:], y[-2:], 'k', alpha=0.1) + self.drawFigure() + + def drawLastPolyLine(self): + self.setXYlim(self.ax.get_xlim(), self.ax.get_ylim()) + x = self._polyx + y = self._polyy + if len(x) >= 2 and len(y) >= 2: + self.ax.plot((x[-1], x[0]), (y[-1], y[0]), 'k', alpha=0.1) + self.drawFigure() + + 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} + self.printOutput('Found %d traces in polygon: %s' % (numtraces, shots)) + + def printOutput(self, text): + print text + self.writeInTextfield(text) + + def chooseRectangles(self, event = None): + ''' + Activates matplotlib widget RectangleSelector. + ''' + from matplotlib.widgets import RectangleSelector + if hasattr(self, '_cidPoly'): + self.disconnectPoly() + self.printOutput('Select rectangle is active. Press and hold left mousebutton.') + self._cidRect = None + self._cidRect = self.ax.figure.canvas.mpl_connect('button_press_event', self._onpress) + self._rectangle = RectangleSelector(self.ax, self._onselect_clicks) + return self._rectangle + + def choosePolygon(self, event = None): + ''' + Activates matplotlib widget LassoSelector. + ''' + from matplotlib.widgets import LassoSelector + if hasattr(self, '_cidRect'): + self.disconnectRect() + self.printOutput('Select polygon is active. Add points with leftclick. Finish with rightclick.') + self._cidPoly = None + self._cidPoly = self.ax.figure.canvas.mpl_connect('button_press_event', self._onpress) + self._lasso = LassoSelector(self.ax, self._onselect_verts) + return self._lasso + + def disconnectPoly(self, event = None): + if not hasattr(self, '_cidPoly'): + self.printOutput('no poly selection found') + return + self.ax.figure.canvas.mpl_disconnect(self._cidPoly) + del self._cidPoly + self.finishPolygon() + self._lasso.disconnect_events() + print 'disconnected poly selection\n' + + def disconnectRect(self, event = None): + if not hasattr(self, '_cidRect'): + self.printOutput('no rectangle selection found') + return + self.ax.figure.canvas.mpl_disconnect(self._cidRect) + del self._cidRect + self._rectangle.disconnect_events() + print 'disconnected rectangle selection\n' + + def deselectLastSelection(self, event = None): + if self.shots_found.keys() == []: + self.printOutput('No selection found.') + return + key = max(self.shots_found.keys()) + self.deselectSelection(key) + + def deselectSelection(self, key, color = 'green', alpha = 0.1): + if key not in self.shots_found.keys(): + self.printOutput('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) + self.printOutput('Deselected selection number %d' % key) + + def findTracesInPoly(self, x, y, picks = 'normal', highlight = True): + def dotproduct(v1, v2): + return sum((a * b for a, b in zip(v1, v2))) + + 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 = 1e-07 + 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: + self.printOutput('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 shotnumber not 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 shotnumber not 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()) + self.ax.annotate(s='s%s|t%s' % (shot.getShotnumber(), traceID), xy=(shot.getDistance(traceID), shot.getPick(traceID)), fontsize='xx-small') - def plotTracesInRegion(self): + def highlightAllActiveRegions(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 plotTracesInActiveRegions(self, event = None, 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: - 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) + for key in keys: + for shotnumber in self.shots_found[key]['shots']: + if shot.getShotnumber() == 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) + self.printOutput('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 setAllActiveRegionsForDeletion(self, event = None): + 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 key in keys: + for shotnumber in self.shots_found[key]['shots'].keys(): + if shotnumber not in self.shots_for_deletion: + self.shots_for_deletion[shotnumber] = [] + for traceID in self.shots_found[key]['shots'][shotnumber]: + if traceID not in self.shots_for_deletion[shotnumber]: + self.shots_for_deletion[shotnumber].append(traceID) + self.deselectSelection(key, color = 'red', alpha = 0.2) - for shotnumber in self.shots_found: - if not shotnumber in self.shots_for_deletion: - self.shots_for_deletion[shotnumber] = [] - for traceID in self.shots_found[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'): + self.printOutput('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 + 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() - for index in range(len(self._getx0())): - x0 = self._getx0()[index] - y0 = self._gety0()[index] - x1 = self._getx1()[index] - y1 = self._gety1()[index] + 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() - self.ax.add_patch(Rectangle((x0, y0), (x1 - x0), (y1 - y0), alpha=0.5, facecolor = color)) - self.refreshFigure() + def clearShotsForDeletion(self): + ''' + Clears the list of shots marked for deletion. + ''' + self.shots_for_deletion = {} + print('Cleared all shots that were set for deletion.') - def markCurrentRegion(self, x0, x1, y0, y1, color = 'grey'): - from matplotlib.patches import Rectangle + def getShotsForDeletion(self): + return self.shots_for_deletion - self.ax.add_patch(Rectangle((x0, y0), (x1 - x0), (y1 - y0), alpha=0.1, facecolor = color)) - self.refreshFigure() + def deleteAllMarkedPicks(self, event = None): + ''' + Deletes all shots set for deletion. + ''' + if len(self.getShotsForDeletion()) is 0: + self.printOutput('No shots set for deletion.') + return - def deleteMarkedPicks(self): 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 setXYlim(self, xlim, ylim): + self._xlim, self._ylim = xlim, ylim + + def refreshLog10SNR(self, event = None): + cbv = 'log10SNR' + self.refreshFigure(self, colorByVal=cbv) + + def refreshPickerror(self, event = None): + cbv = 'pickerror' + self.refreshFigure(self, colorByVal=cbv) + + def refreshSPE(self, event = None): + cbv = 'spe' + self.refreshFigure(self, colorByVal=cbv) + + def refreshFigure(self, event = None, colorByVal = None): + if colorByVal == None: + colorByVal = self.cbv + else: + self.cbv = colorByVal + self.printOutput('Refreshing figure...') + self.ax.clear() + self.ax = self.survey.plotAllPicks(ax=self.ax, cbar=self.cbar, refreshPlot=True, colorByVal=colorByVal) + self.setXYlim(self.ax.get_xlim(), self.ax.get_ylim()) + self.markAllActiveRegions() + self.highlightAllActiveRegions() + self.drawFigure() + self.printOutput('Done!') + + def drawFigure(self, resetAxes = True): + if resetAxes == True: + self.ax.set_xlim(self._xlim) + self.ax.set_ylim(self._ylim) plt.draw() diff --git a/pylot/core/active/surveyUtils.py b/pylot/core/active/surveyUtils.py index 909fed13..d6fb1d15 100644 --- a/pylot/core/active/surveyUtils.py +++ b/pylot/core/active/surveyUtils.py @@ -1,19 +1,21 @@ -#!/usr/bin/env python -# -*- coding: utf-8 -*- - import numpy as np def readParameters(parfile, parameter): from ConfigParser import ConfigParser parameterConfig = ConfigParser() parameterConfig.read('parfile') - - value = parameterConfig.get('vars', parameter).split('#')[0] - value = value.replace(" ", "") + + value = parameterConfig.get('vars', parameter).split('\t')[0] return value +def setArtificialPick(shot_dict, traceID, pick): + for shot in shot_dict.values(): + shot.setPick(traceID, pick) + shot.setPickwindow(traceID, shot.getCut()) + def fitSNR4dist(shot_dict, shiftdist = 5): + import numpy as np dists = [] picks = [] snrs = [] @@ -23,7 +25,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) @@ -34,6 +36,7 @@ def fitSNR4dist(shot_dict, shiftdist = 5): plotFittedSNR(dists, snrthresholds, snrs) return fit_fn #### ZU VERBESSERN, sollte fertige funktion wiedergeben + def plotFittedSNR(dists, snrthresholds, snrs): import matplotlib.pyplot as plt plt.interactive(True) @@ -45,25 +48,28 @@ def plotFittedSNR(dists, snrthresholds, snrs): plt.legend() def setFittedSNR(shot_dict, shiftdist = 5, p1 = 0.004, p2 = -0.004): + import numpy as np #fit_fn = fitSNR4dist(shot_dict) fit_fn = np.poly1d([p1, p2]) for shot in shot_dict.values(): for traceID in shot.getTraceIDlist(): ### IMPROVE shot.setSNRthreshold(traceID, 1/(fit_fn(shot.getDistance(traceID) + shiftdist)**2)) ### s.o. - print "\nsetFittedSNR: Finished setting of fitted SNR-threshold" + print "setFittedSNR: Finished setting of fitted SNR-threshold" + def findTracesInRanges(shot_dict, distancebin, pickbin): ''' Returns traces corresponding to a certain area in a plot with all picks over the distances. - + :param: shot_dict, dictionary containing all shots that are used :type: dictionary - + :param: distancebin :type: tuple, (dist1[m], dist2[m]) - + :param: pickbin :type: tuple, (t1[s], t2[s]) + ''' shots_found = {} for shot in shot_dict.values(): diff --git a/pylot/core/analysis/__init__.py b/pylot/core/analysis/__init__.py index e69de29b..40a96afc 100644 --- a/pylot/core/analysis/__init__.py +++ b/pylot/core/analysis/__init__.py @@ -0,0 +1 @@ +# -*- coding: utf-8 -*- diff --git a/pylot/core/analysis/coinctimes.py b/pylot/core/analysis/coinctimes.py index b1486898..1ecc7a9e 100644 --- a/pylot/core/analysis/coinctimes.py +++ b/pylot/core/analysis/coinctimes.py @@ -6,7 +6,7 @@ from obspy.signal.trigger import coincidenceTrigger -class CoincidenceTimes(): +class CoincidenceTimes(object): def __init__(self, st, comp='Z', coinum=4, sta=1., lta=10., on=5., off=1.): _type = 'recstalta' @@ -54,4 +54,4 @@ def main(): if __name__ == '__main__': - main() \ No newline at end of file + main() diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index ab49bf86..ab455762 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -1,3 +1,4 @@ +#!/usr/bin/env python # -*- coding: utf-8 -*- """ Created August/September 2015. @@ -34,7 +35,7 @@ class Magnitude(object): :type: integer ''' - + assert isinstance(wfstream, Stream), "%s is not a stream object" % str(wfstream) self.setwfstream(wfstream) @@ -62,7 +63,7 @@ class Magnitude(object): def setpwin(self, pwin): self.pwin = pwin - + def getiplot(self): return self.iplot @@ -71,7 +72,7 @@ class Magnitude(object): def getwapp(self): return self.wapp - + def getw0(self): return self.w0 @@ -103,7 +104,7 @@ class WApp(Magnitude): 'poles': [5.6089 - 5.4978j, -5.6089 - 5.4978j], 'zeros': [0j, 0j], 'gain': 2080, - 'sensitivity': 1} + 'sensitivity': 1} stream.simulate(paz_remove=None, paz_simulate=paz_wa) @@ -133,19 +134,19 @@ class WApp(Magnitude): raw_input() plt.close(f) - + class DCfc(Magnitude): ''' - Method to calculate the source spectrum and to derive from that the plateau - (so-called DC-value) and the corner frequency assuming Aki's omega-square + Method to calculate the source spectrum and to derive from that the plateau + (so-called DC-value) and the corner frequency assuming Aki's omega-square source model. Has to be derived from instrument corrected displacement traces! ''' def calcsourcespec(self): - print ("Calculating source spectrum ....") + print ("Calculating source spectrum ....") self.w0 = None # DC-value - self.fc = None # corner frequency + self.fc = None # corner frequency stream = self.getwfstream() tr = stream[0] @@ -155,7 +156,7 @@ class DCfc(Magnitude): iwin = getsignalwin(t, self.getTo(), self.getpwin()) xdat = tr.data[iwin] - # fft + # fft fny = tr.stats.sampling_rate / 2 l = len(xdat) / tr.stats.sampling_rate n = tr.stats.sampling_rate * l # number of fft bins after Bath @@ -167,7 +168,7 @@ class DCfc(Magnitude): L = (N - 1) / tr.stats.sampling_rate f = np.arange(0, fny, 1/L) - # remove zero-frequency and frequencies above + # remove zero-frequency and frequencies above # corner frequency of seismometer (assumed # to be 100 Hz) fi = np.where((f >= 1) & (f < 100)) @@ -184,18 +185,16 @@ class DCfc(Magnitude): [optspecfit, pcov] = curve_fit(synthsourcespec, F, YY.real, [DCin, Fcin]) self.w0 = optspecfit[0] self.fc = optspecfit[1] - print ("DCfc: Determined DC-value: %e m/Hz, \n" \ - "Determined corner frequency: %f Hz" % (self.w0, self.fc)) - - - #if self.getiplot() > 1: - iplot=2 - if iplot > 1: + print ("DCfc: Determined DC-value: %e m/Hz, \n" + "Determined corner frequency: %f Hz" % (self.w0, self.fc)) + + + if self.getiplot() > 1: f1 = plt.figure() plt.subplot(2,1,1) # show displacement in mm - plt.plot(t, np.multiply(tr, 1000), 'k') - plt.plot(t[iwin], np.multiply(xdat, 1000), 'g') + plt.plot(t, np.multiply(tr, 1000), 'k') + plt.plot(t[iwin], np.multiply(xdat, 1000), 'g') plt.title('Seismogram and P pulse, station %s' % tr.stats.station) plt.xlabel('Time since %s' % tr.stats.starttime) plt.ylabel('Displacement [mm]') @@ -216,8 +215,8 @@ class DCfc(Magnitude): def synthsourcespec(f, omega0, fcorner): ''' - Calculates synthetic source spectrum from given plateau and corner - frequency assuming Akis omega-square model. + Calculates synthetic source spectrum from given plateau and corner + frequency assuming Akis omega-square model. :param: f, frequencies :type: array @@ -228,7 +227,7 @@ def synthsourcespec(f, omega0, fcorner): :param: fcorner, corner frequency of source spectrum :type: float ''' - + #ssp = omega0 / (pow(2, (1 + f / fcorner))) ssp = omega0 / (1 + pow(2, (f / fcorner))) diff --git a/pylot/core/pick/CharFuns.py b/pylot/core/pick/CharFuns.py index 2f6e8cd8..78efd083 100644 --- a/pylot/core/pick/CharFuns.py +++ b/pylot/core/pick/CharFuns.py @@ -1,3 +1,4 @@ +#!/usr/bin/env python # -*- coding: utf-8 -*- """ Created Oct/Nov 2014 @@ -24,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. @@ -61,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 @@ -119,7 +121,7 @@ class CharacteristicFunction(object): def getTimeArray(self): incr = self.getIncrement() - self.TimeArray = np.arange(0, len(self.getCF()) * incr, incr) + self.getCut()[0] + self.TimeArray = np.arange(0, len(self.getCF()) * incr, incr) + self.getCut()[0] return self.TimeArray def getFnoise(self): @@ -134,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) @@ -164,31 +169,31 @@ class CharacteristicFunction(object): stop = min([len(self.orig_data[0]), len(self.orig_data[1])]) elif self.cut[0] == 0 and self.cut[1] is not 0: start = 0 - stop = min([self.cut[1] / self.dt, len(self.orig_data[0]), \ - len(self.orig_data[1])]) + stop = min([self.cut[1] / self.dt, len(self.orig_data[0]), + len(self.orig_data[1])]) else: start = max([0, self.cut[0] / self.dt]) - stop = min([self.cut[1] / self.dt, len(self.orig_data[0]), \ - len(self.orig_data[1])]) + stop = min([self.cut[1] / self.dt, len(self.orig_data[0]), + len(self.orig_data[1])]) hh = self.orig_data.copy() h1 = hh[0].copy() h2 = hh[1].copy() hh[0].data = h1.data[int(start):int(stop)] hh[1].data = h2.data[int(start):int(stop)] - data = hh + data = hh return data elif len(self.orig_data) == 3: if self.cut[0] == 0 and self.cut[1] == 0: start = 0 - stop = min([self.cut[1] / self.dt, len(self.orig_data[0]), \ - len(self.orig_data[1]), len(self.orig_data[2])]) + stop = min([self.cut[1] / self.dt, len(self.orig_data[0]), + len(self.orig_data[1]), len(self.orig_data[2])]) elif self.cut[0] == 0 and self.cut[1] is not 0: start = 0 stop = self.cut[1] / self.dt else: start = max([0, self.cut[0] / self.dt]) - stop = min([self.cut[1] / self.dt, len(self.orig_data[0]), \ - len(self.orig_data[1]), len(self.orig_data[2])]) + stop = min([self.cut[1] / self.dt, len(self.orig_data[0]), + len(self.orig_data[1]), len(self.orig_data[2])]) hh = self.orig_data.copy() h1 = hh[0].copy() h2 = hh[1].copy() @@ -196,12 +201,12 @@ class CharacteristicFunction(object): hh[0].data = h1.data[int(start):int(stop)] hh[1].data = h2.data[int(start):int(stop)] hh[2].data = h3.data[int(start):int(stop)] - data = hh + data = hh return data else: data = self.orig_data.copy() return data - + def calcCF(self, data=None): self.cf = data @@ -218,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) @@ -230,7 +236,7 @@ class AICcf(CharacteristicFunction): cumsumcf = np.cumsum(np.power(xnp, 2)) i = np.where(cumsumcf == 0) cumsumcf[i] = np.finfo(np.float64).eps - cf[k] = ((k - 1) * np.log(cumsumcf[k] / k) + (datlen - k + 1) * \ + cf[k] = ((k - 1) * np.log(cumsumcf[k] / k) + (datlen - k + 1) * np.log((cumsumcf[datlen - 1] - cumsumcf[k - 1]) / (datlen - k + 1))) cf[0] = cf[1] inf = np.isinf(cf) @@ -256,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) @@ -285,7 +293,7 @@ class HOScf(CharacteristicFunction): LTA[j] = lta / np.power(lta1, 1.5) elif self.getOrder() == 4: LTA[j] = lta / np.power(lta1, 2) - + nn = np.isnan(LTA) if len(nn) > 1: LTA[nn] = 0 @@ -315,7 +323,7 @@ class ARZcf(CharacteristicFunction): cf = np.zeros(len(xnp)) loopstep = self.getARdetStep() arcalci = ldet + self.getOrder() #AR-calculation index - for i in range(ldet + self.getOrder(), tend - lpred - 1): + for i in range(ldet + self.getOrder(), tend - lpred - 1): if i == arcalci: #determination of AR coefficients #to speed up calculation, AR-coefficients are calculated only every i+loopstep[1]! @@ -362,7 +370,7 @@ class ARZcf(CharacteristicFunction): rhs = np.zeros(self.getOrder()) for k in range(0, self.getOrder()): for i in range(rind, ldet+1): - ki = k + 1 + ki = k + 1 rhs[k] = rhs[k] + data[i] * data[i - ki] #recursive calculation of data array (second sum at left part of eq. 6.5 in Kueperkoch et al. 2012) @@ -382,7 +390,7 @@ class ARZcf(CharacteristicFunction): def arPredZ(self, data, arpara, rind, lpred): ''' Function to predict waveform, assuming an autoregressive process of order - p (=size(arpara)), with AR parameters arpara calculated in arDet. After + p (=size(arpara)), with AR parameters arpara calculated in arDet. After Thomas Meier (CAU), published in Kueperkoch et al. (2012). :param: data, time series to be predicted :type: array @@ -400,9 +408,9 @@ class ARZcf(CharacteristicFunction): ''' #be sure of the summation indeces if rind < len(arpara): - rind = len(arpara) + rind = len(arpara) if rind > len(data) - lpred : - rind = len(data) - lpred + rind = len(data) - lpred if lpred < 1: lpred = 1 if lpred > len(data) - 2: @@ -422,7 +430,7 @@ class ARHcf(CharacteristicFunction): def calcCF(self, data): print 'Calculating AR-prediction error from both horizontal traces ...' - + xnp = self.getDataArray(self.getCut()) n0 = np.isnan(xnp[0].data) if len(n0) > 1: @@ -430,7 +438,7 @@ class ARHcf(CharacteristicFunction): n1 = np.isnan(xnp[1].data) if len(n1) > 1: xnp[1].data[n1] = 0 - + #some parameters needed #add noise to time series xenoise = xnp[0].data + np.random.normal(0.0, 1.0, len(xnp[0].data)) * self.getFnoise() * max(abs(xnp[0].data)) @@ -441,7 +449,7 @@ class ARHcf(CharacteristicFunction): #Time2: length of AR-prediction window [sec] ldet = int(round(self.getTime1() / self.getIncrement())) #length of AR-determination window [samples] lpred = int(np.ceil(self.getTime2() / self.getIncrement())) #length of AR-prediction window [samples] - + cf = np.zeros(len(xenoise)) loopstep = self.getARdetStep() arcalci = lpred + self.getOrder() - 1 #AR-calculation index @@ -515,7 +523,7 @@ class ARHcf(CharacteristicFunction): def arPredH(self, data, arpara, rind, lpred): ''' Function to predict waveform, assuming an autoregressive process of order - p (=size(arpara)), with AR parameters arpara calculated in arDet. After + p (=size(arpara)), with AR parameters arpara calculated in arDet. After Thomas Meier (CAU), published in Kueperkoch et al. (2012). :param: data, horizontal component seismograms to be predicted :type: structured array @@ -558,7 +566,7 @@ class AR3Ccf(CharacteristicFunction): def calcCF(self, data): print 'Calculating AR-prediction error from all 3 components ...' - + xnp = self.getDataArray(self.getCut()) n0 = np.isnan(xnp[0].data) if len(n0) > 1: @@ -569,7 +577,7 @@ class AR3Ccf(CharacteristicFunction): n2 = np.isnan(xnp[2].data) if len(n2) > 1: xnp[2].data[n2] = 0 - + #some parameters needed #add noise to time series xenoise = xnp[0].data + np.random.normal(0.0, 1.0, len(xnp[0].data)) * self.getFnoise() * max(abs(xnp[0].data)) @@ -581,7 +589,7 @@ class AR3Ccf(CharacteristicFunction): #Time2: length of AR-prediction window [sec] ldet = int(round(self.getTime1() / self.getIncrement())) #length of AR-determination window [samples] lpred = int(np.ceil(self.getTime2() / self.getIncrement())) #length of AR-prediction window [samples] - + cf = np.zeros(len(xenoise)) loopstep = self.getARdetStep() arcalci = ldet + self.getOrder() - 1 #AR-calculation index @@ -616,7 +624,7 @@ class AR3Ccf(CharacteristicFunction): Function to calculate AR parameters arpara after Thomas Meier (CAU), published in Kueperkoch et al. (2012). This function solves SLE using the Moore- Penrose inverse, i.e. the least-squares approach. "data" is a structured array. - AR parameters are calculated based on both horizontal components and vertical + AR parameters are calculated based on both horizontal components and vertical componant. :param: data, horizontal component seismograms to calculate AR parameters from :type: structured array @@ -658,7 +666,7 @@ class AR3Ccf(CharacteristicFunction): def arPred3C(self, data, arpara, rind, lpred): ''' Function to predict waveform, assuming an autoregressive process of order - p (=size(arpara)), with AR parameters arpara calculated in arDet3C. After + p (=size(arpara)), with AR parameters arpara calculated in arDet3C. After Thomas Meier (CAU), published in Kueperkoch et al. (2012). :param: data, horizontal and vertical component seismograms to be predicted :type: structured array diff --git a/pylot/core/pick/Picker.py b/pylot/core/pick/Picker.py index f9ef82d4..41ac2439 100644 --- a/pylot/core/pick/Picker.py +++ b/pylot/core/pick/Picker.py @@ -312,7 +312,7 @@ class PragPicker(AutoPicking): else: for i in range(1, len(self.cf)): if i > ismooth: - ii1 = i - ismooth; + ii1 = i - ismooth cfsmooth[i] = cfsmooth[i - 1] + (self.cf[i] - self.cf[ii1]) / ismooth else: cfsmooth[i] = np.mean(self.cf[1 : i]) diff --git a/pylot/core/pick/__init__.py b/pylot/core/pick/__init__.py index 4287ca86..ec51c5a2 100644 --- a/pylot/core/pick/__init__.py +++ b/pylot/core/pick/__init__.py @@ -1 +1,2 @@ -# \ No newline at end of file +# -*- coding: utf-8 -*- +# diff --git a/pylot/core/pick/autopick.py b/pylot/core/pick/autopick.py index e3771d7a..19d33b02 100755 --- a/pylot/core/pick/autopick.py +++ b/pylot/core/pick/autopick.py @@ -204,27 +204,27 @@ def autopickstation(wfstream, pickparam): if len(ndat) == 0 or len(edat) == 0: print ("One or more horizontal components missing!") print ("Signal length only checked on vertical component!") - print ("Decreasing minsiglengh from %f to %f" \ - % (minsiglength, minsiglength / 2)) + print ("Decreasing minsiglengh from %f to %f" + % (minsiglength, minsiglength / 2)) Pflag = checksignallength(zne, aicpick.getpick(), tsnrz, - minsiglength / 2, \ + minsiglength / 2, nfacsl, minpercent, iplot) else: # filter and taper horizontal traces trH1_filt = edat.copy() trH2_filt = ndat.copy() trH1_filt.filter('bandpass', freqmin=bph1[0], - freqmax=bph1[1], \ - zerophase=False) + freqmax=bph1[1], + zerophase=False) trH2_filt.filter('bandpass', freqmin=bph1[0], - freqmax=bph1[1], \ - zerophase=False) + freqmax=bph1[1], + zerophase=False) trH1_filt.taper(max_percentage=0.05, type='hann') trH2_filt.taper(max_percentage=0.05, type='hann') zne += trH1_filt zne += trH2_filt Pflag = checksignallength(zne, aicpick.getpick(), tsnrz, - minsiglength, \ + minsiglength, nfacsl, minpercent, iplot) if Pflag == 1: @@ -234,7 +234,7 @@ def autopickstation(wfstream, pickparam): print 'One or more horizontal components missing!' print 'Skipping control function checkZ4S.' else: - Pflag = checkZ4S(zne, aicpick.getpick(), zfac, \ + Pflag = checkZ4S(zne, aicpick.getpick(), zfac, tsnrz[3], iplot) if Pflag == 0: Pmarker = 'SinsteadP' @@ -317,31 +317,31 @@ def autopickstation(wfstream, pickparam): data = Data() [corzdat, restflag] = data.restituteWFData(invdir, zdat) if restflag == 1: - # integrate to displacement - corintzdat = integrate.cumtrapz(corzdat[0], None, corzdat[0].stats.delta) - # class needs stream object => build it - z_copy = zdat.copy() - z_copy[0].data = corintzdat - # largest detectable period == window length - # after P pulse for calculating source spectrum - winzc = (1 / bpz2[0]) * z_copy[0].stats.sampling_rate - impickP = mpickP * z_copy[0].stats.sampling_rate - wfzc = z_copy[0].data[impickP : impickP + winzc] - # calculate spectrum using only first cycles of - # waveform after P onset! - zc = crossings_nonzero_all(wfzc) - if np.size(zc) == 0: - print ("Something is wrong with the waveform, " \ - "no zero crossings derived!") - print ("Cannot calculate source spectrum!") - else: - calcwin = (zc[3] - zc[0]) * z_copy[0].stats.delta - # calculate source spectrum and get w0 and fc - specpara = DCfc(z_copy, mpickP, calcwin, iplot) - w0 = specpara.getw0() - fc = specpara.getfc() + # integrate to displacement + corintzdat = integrate.cumtrapz(corzdat[0], None, corzdat[0].stats.delta) + # class needs stream object => build it + z_copy = zdat.copy() + z_copy[0].data = corintzdat + # largest detectable period == window length + # after P pulse for calculating source spectrum + winzc = (1 / bpz2[0]) * z_copy[0].stats.sampling_rate + impickP = mpickP * z_copy[0].stats.sampling_rate + wfzc = z_copy[0].data[impickP : impickP + winzc] + # calculate spectrum using only first cycles of + # waveform after P onset! + zc = crossings_nonzero_all(wfzc) + if np.size(zc) == 0: + print ("Something is wrong with the waveform, " + "no zero crossings derived!") + print ("Cannot calculate source spectrum!") + else: + calcwin = (zc[3] - zc[0]) * z_copy[0].stats.delta + # calculate source spectrum and get w0 and fc + specpara = DCfc(z_copy, mpickP, calcwin, iplot) + w0 = specpara.getw0() + fc = specpara.getfc() - print ("autopickstation: P-weight: %d, SNR: %f, SNR[dB]: %f, " \ + print ("autopickstation: P-weight: %d, SNR: %f, SNR[dB]: %f, " "Polarity: %s" % (Pweight, SNRP, SNRPdB, FM)) Sflag = 1 @@ -352,7 +352,7 @@ def autopickstation(wfstream, pickparam): Sflag = 0 else: - print ("autopickstation: No vertical component data available!, " \ + print ("autopickstation: No vertical component data available!, " "Skipping station!") if edat is not None and ndat is not None and len(edat) > 0 and len( @@ -560,7 +560,7 @@ def autopickstation(wfstream, pickparam): hdat += ndat h_copy = hdat.copy() [cordat, restflag] = data.restituteWFData(invdir, h_copy) - # calculate WA-peak-to-peak amplitude + # calculate WA-peak-to-peak amplitude # using subclass WApp of superclass Magnitude if restflag == 1: if Sweight < 4: @@ -591,10 +591,10 @@ def autopickstation(wfstream, pickparam): h_copy = hdat.copy() [cordat, restflag] = data.restituteWFData(invdir, h_copy) if restflag == 1: - # calculate WA-peak-to-peak amplitude + # calculate WA-peak-to-peak amplitude # using subclass WApp of superclass Magnitude - wapp = WApp(cordat, mpickP, mpickP + sstop + (0.5 * (mpickP \ - + sstop)), iplot) + wapp = WApp(cordat, mpickP, mpickP + sstop + (0.5 * (mpickP + + sstop)), iplot) Ao = wapp.getwapp() else: @@ -771,14 +771,14 @@ def autopickstation(wfstream, pickparam): # create dictionary # for P phase phase = 'P' - phasepick = {'lpp': lpickP, 'epp': epickP, 'mpp': mpickP, 'spe': Perror, \ + phasepick = {'lpp': lpickP, 'epp': epickP, 'mpp': mpickP, 'spe': Perror, 'snr': SNRP, 'snrdb': SNRPdB, 'weight': Pweight, 'fm': FM} picks = {phase: phasepick} # add P marker picks[phase]['marked'] = Pmarker # add S phase phase = 'S' - phasepick = {'lpp': lpickS, 'epp': epickS, 'mpp': mpickS, 'spe': Serror, \ + phasepick = {'lpp': lpickS, 'epp': epickS, 'mpp': mpickS, 'spe': Serror, 'snr': SNRS, 'snrdb': SNRSdB, 'weight': Sweight, 'fm': None} picks[phase] = phasepick # add Wood-Anderson amplitude diff --git a/pylot/core/pick/run_makeCF.py b/pylot/core/pick/run_makeCF.py index d1748d67..2ecd636f 100755 --- a/pylot/core/pick/run_makeCF.py +++ b/pylot/core/pick/run_makeCF.py @@ -6,8 +6,8 @@ Only for test purposes! """ -from obspy.core import read -import matplotlib.pyplot as plt +from obspy.core import read +import matplotlib.pyplot as plt import numpy as np from pylot.core.pick.CharFuns import CharacteristicFunction from pylot.core.pick.Picker import AutoPicking @@ -56,7 +56,7 @@ def run_makeCF(project, database, event, iplot, station=None): st_copy = st.copy() #filter and taper data tr_filt = st[0].copy() - tr_filt.filter('bandpass', freqmin=bpz[0], freqmax=bpz[1], zerophase=False) + tr_filt.filter('bandpass', freqmin=bpz[0], freqmax=bpz[1], zerophase=False) tr_filt.taper(max_percentage=0.05, type='hann') st_copy[0].data = tr_filt.data ############################################################## @@ -120,8 +120,8 @@ def run_makeCF(project, database, event, iplot, station=None): #filter and taper data trH1_filt = H[0].copy() trH2_filt = H[1].copy() - trH1_filt.filter('bandpass', freqmin=bph[0], freqmax=bph[1], zerophase=False) - trH2_filt.filter('bandpass', freqmin=bph[0], freqmax=bph[1], zerophase=False) + trH1_filt.filter('bandpass', freqmin=bph[0], freqmax=bph[1], zerophase=False) + trH2_filt.filter('bandpass', freqmin=bph[0], freqmax=bph[1], zerophase=False) trH1_filt.taper(max_percentage=0.05, type='hann') trH2_filt.taper(max_percentage=0.05, type='hann') H_copy[0].data = trH1_filt.data @@ -167,9 +167,9 @@ def run_makeCF(project, database, event, iplot, station=None): All1_filt = AllC[0].copy() All2_filt = AllC[1].copy() All3_filt = AllC[2].copy() - All1_filt.filter('bandpass', freqmin=bph[0], freqmax=bph[1], zerophase=False) - All2_filt.filter('bandpass', freqmin=bph[0], freqmax=bph[1], zerophase=False) - All3_filt.filter('bandpass', freqmin=bpz[0], freqmax=bpz[1], zerophase=False) + All1_filt.filter('bandpass', freqmin=bph[0], freqmax=bph[1], zerophase=False) + All2_filt.filter('bandpass', freqmin=bph[0], freqmax=bph[1], zerophase=False) + All3_filt.filter('bandpass', freqmin=bpz[0], freqmax=bpz[1], zerophase=False) All1_filt.taper(max_percentage=0.05, type='hann') All2_filt.taper(max_percentage=0.05, type='hann') All3_filt.taper(max_percentage=0.05, type='hann') @@ -209,19 +209,19 @@ def run_makeCF(project, database, event, iplot, station=None): plt.ylim([-1.5, 1.5]) plt.xlabel('Time [s]') plt.ylabel('Normalized Counts') - plt.title('%s, %s, CF-SNR=%7.2f, CF-Slope=%12.2f' % (tr.stats.station, \ - tr.stats.channel, aicpick.getSNR(), aicpick.getSlope())) + plt.title('%s, %s, CF-SNR=%7.2f, CF-Slope=%12.2f' % (tr.stats.station, + tr.stats.channel, aicpick.getSNR(), aicpick.getSlope())) plt.suptitle(tr.stats.starttime) - plt.legend([p1, p2, p3, p4, p5], ['Data', 'HOS-CF', 'HOSAIC-CF', 'ARZ-CF', 'ARZAIC-CF']) + plt.legend([p1, p2, p3, p4, p5], ['Data', 'HOS-CF', 'HOSAIC-CF', 'ARZ-CF', 'ARZAIC-CF']) #plot horizontal traces plt.figure(2) plt.subplot(2,1,1) - tsteph = tpredh / 4 + tsteph = tpredh / 4 th1data = np.arange(0, trH1_filt.stats.npts / trH1_filt.stats.sampling_rate, trH1_filt.stats.delta) th2data = np.arange(0, trH2_filt.stats.npts / trH2_filt.stats.sampling_rate, trH2_filt.stats.delta) tarhcf = np.arange(0, len(arhcf.getCF()) * tsteph, tsteph) + cuttimes[0] + tdeth +tpredh p21, = plt.plot(th1data, trH1_filt.data/max(trH1_filt.data), 'k') - p22, = plt.plot(arhcf.getTimeArray(), arhcf.getCF()/max(arhcf.getCF()), 'r') + p22, = plt.plot(arhcf.getTimeArray(), arhcf.getCF()/max(arhcf.getCF()), 'r') p23, = plt.plot(arhaiccf.getTimeArray(), arhaiccf.getCF()/max(arhaiccf.getCF())) plt.plot([aicarhpick.getpick(), aicarhpick.getpick()], [-1, 1], 'b') plt.plot([aicarhpick.getpick()-0.5, aicarhpick.getpick()+0.5], [1, 1], 'b') @@ -238,10 +238,10 @@ def run_makeCF(project, database, event, iplot, station=None): plt.ylabel('Normalized Counts') plt.title([trH1_filt.stats.station, trH1_filt.stats.channel]) plt.suptitle(trH1_filt.stats.starttime) - plt.legend([p21, p22, p23], ['Data', 'ARH-CF', 'ARHAIC-CF']) + plt.legend([p21, p22, p23], ['Data', 'ARH-CF', 'ARHAIC-CF']) plt.subplot(2,1,2) plt.plot(th2data, trH2_filt.data/max(trH2_filt.data), 'k') - plt.plot(arhcf.getTimeArray(), arhcf.getCF()/max(arhcf.getCF()), 'r') + plt.plot(arhcf.getTimeArray(), arhcf.getCF()/max(arhcf.getCF()), 'r') plt.plot(arhaiccf.getTimeArray(), arhaiccf.getCF()/max(arhaiccf.getCF())) plt.plot([aicarhpick.getpick(), aicarhpick.getpick()], [-1, 1], 'b') plt.plot([aicarhpick.getpick()-0.5, aicarhpick.getpick()+0.5], [1, 1], 'b') @@ -271,7 +271,7 @@ def run_makeCF(project, database, event, iplot, station=None): plt.ylabel('Normalized Counts') plt.title([tr.stats.station, tr.stats.channel]) plt.suptitle(trH1_filt.stats.starttime) - plt.legend([p31, p32], ['Data', 'AR3C-CF']) + plt.legend([p31, p32], ['Data', 'AR3C-CF']) plt.subplot(3,1,2) plt.plot(th1data, trH1_filt.data/max(trH1_filt.data), 'k') plt.plot(ar3ccf.getTimeArray(), ar3ccf.getCF()/max(ar3ccf.getCF()), 'r') diff --git a/pylot/core/pick/utils.py b/pylot/core/pick/utils.py index 233bb70d..9feabbb8 100644 --- a/pylot/core/pick/utils.py +++ b/pylot/core/pick/utils.py @@ -1,4 +1,5 @@ #!/usr/bin/env python +# -*- coding: utf-8 -*- # # -*- coding: utf-8 -*- """ @@ -14,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, @@ -43,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, @@ -74,8 +76,9 @@ 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: - print("earllatepicker: Doubled signal window size %s time(s) " - "because of NaN for earliest pick." %count) + if stealthMode is False: + print("\nearllatepicker: Doubled signal window size %s time(s) " + "because of NaN for earliest pick." %count) isigDoubleWinStart = pis[-1] + 1 isignalDoubleWin = np.arange(isigDoubleWinStart, isigDoubleWinStart + len(pis)) @@ -91,7 +94,7 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None): T0 = np.mean(np.diff(zc)) * X[0].stats.delta # this is half wave length # T0/4 is assumed as time difference between most likely and earliest possible pick! EPick = Pick1 - T0 / 2 - + # get symmetric pick error as mean from earliest and latest possible pick # by weighting latest possible pick two times earliest possible pick @@ -109,7 +112,7 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None): markersize=14) plt.legend([p1, p2, p3, p4, p5], ['Data', 'Noise Window', 'Signal Window', 'Noise Level', - 'Zero Crossings'], \ + 'Zero Crossings'], loc='best') plt.plot([t[0], t[int(len(t)) - 1]], [-nlevel, -nlevel], '--k') plt.plot([Pick1, Pick1], [max(x), -max(x)], 'b', linewidth=2) @@ -182,10 +185,10 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None): i = 0 for j in range(ipick[0][1], ipick[0][len(t[ipick]) - 1]): i = i + 1 - if xraw[j - 1] <= 0 and xraw[j] >= 0: + if xraw[j - 1] <= 0 <= xraw[j]: zc1.append(t[ipick][i]) index1.append(i) - elif xraw[j - 1] > 0 and xraw[j] <= 0: + elif xraw[j - 1] > 0 >= xraw[j]: zc1.append(t[ipick][i]) index1.append(i) if len(zc1) == 3: @@ -224,10 +227,10 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None): i = 0 for j in range(ipick[0][1], ipick[0][len(t[ipick]) - 1]): i = i + 1 - if xfilt[j - 1] <= 0 and xfilt[j] >= 0: + if xfilt[j - 1] <= 0 <= xfilt[j]: zc2.append(t[ipick][i]) index2.append(i) - elif xfilt[j - 1] > 0 and xfilt[j] <= 0: + elif xfilt[j - 1] > 0 >= xfilt[j]: zc2.append(t[ipick][i]) index2.append(i) if len(zc2) == 3: @@ -262,15 +265,15 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None): if P1 is not None and P2 is not None: if P1[0] < 0 and P2[0] < 0: FM = 'D' - elif P1[0] >= 0 and P2[0] < 0: + elif P1[0] >= 0 > P2[0]: FM = '-' - elif P1[0] < 0 and P2[0] >= 0: + elif P1[0] < 0 <= P2[0]: FM = '-' elif P1[0] > 0 and P2[0] > 0: FM = 'U' - elif P1[0] <= 0 and P2[0] > 0: + elif P1[0] <= 0 < P2[0]: FM = '+' - elif P1[0] > 0 and P2[0] <= 0: + elif P1[0] > 0 >= P2[0]: FM = '+' print ("fmpicker: Found polarity %s" % FM) @@ -285,7 +288,7 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None): p3, = plt.plot(zc1, np.zeros(len(zc1)), '*g', markersize=14) p4, = plt.plot(t[islope1], datafit1, '--g', linewidth=2) plt.legend([p1, p2, p3, p4], - ['Pick', 'Slope Window', 'Zero Crossings', 'Slope'], \ + ['Pick', 'Slope Window', 'Zero Crossings', 'Slope'], loc='best') plt.text(Pick + 0.02, max(xraw) / 2, '%s' % FM, fontsize=14) ax = plt.gca() @@ -493,9 +496,9 @@ def wadaticheck(pickdic, dttolerance, iplot): if len(SPtimes) >= 3: - # calculate slope - p1 = np.polyfit(Ppicks, SPtimes, 1) - wdfit = np.polyval(p1, Ppicks) + # calculate slope + p1 = np.polyfit(Ppicks, SPtimes, 1) + wdfit = np.polyval(p1, Ppicks) wfitflag = 0 # calculate vp/vs ratio before check @@ -532,40 +535,40 @@ def wadaticheck(pickdic, dttolerance, iplot): pickdic[key]['S']['marked'] = marker if len(checkedPpicks) >= 3: - # calculate new slope - p2 = np.polyfit(checkedPpicks, checkedSPtimes, 1) - wdfit2 = np.polyval(p2, checkedPpicks) + # calculate new slope + p2 = np.polyfit(checkedPpicks, checkedSPtimes, 1) + wdfit2 = np.polyval(p2, checkedPpicks) - # calculate vp/vs ratio after check - cvpvsr = p2[0] + 1 - print ("wadaticheck: Average Vp/Vs ratio after check: %f" % cvpvsr) - print ("wadatacheck: Skipped %d S pick(s)" % ibad) + # calculate vp/vs ratio after check + cvpvsr = p2[0] + 1 + print ("wadaticheck: Average Vp/Vs ratio after check: %f" % cvpvsr) + print ("wadatacheck: Skipped %d S pick(s)" % ibad) else: - print ("###############################################") - print ("wadatacheck: Not enough checked S-P times available!") - print ("Skip Wadati check!") + print ("###############################################") + print ("wadatacheck: Not enough checked S-P times available!") + print ("Skip Wadati check!") checkedonsets = pickdic else: - print ("wadaticheck: Not enough S-P times available for reliable regression!") + print ("wadaticheck: Not enough S-P times available for reliable regression!") print ("Skip wadati check!") wfitflag = 1 # plot results if iplot > 1: - plt.figure(iplot) - f1, = plt.plot(Ppicks, SPtimes, 'ro') + plt.figure(iplot) + f1, = plt.plot(Ppicks, SPtimes, 'ro') if wfitflag == 0: - f2, = plt.plot(Ppicks, wdfit, 'k') - f3, = plt.plot(checkedPpicks, checkedSPtimes, 'ko') - f4, = plt.plot(checkedPpicks, wdfit2, 'g') - plt.title('Wadati-Diagram, %d S-P Times, Vp/Vs(raw)=%5.2f,' \ - 'Vp/Vs(checked)=%5.2f' % (len(SPtimes), vpvsr, cvpvsr)) - plt.legend([f1, f2, f3, f4], ['Skipped S-Picks', 'Wadati 1', \ - 'Reliable S-Picks', 'Wadati 2'], loc='best') + f2, = plt.plot(Ppicks, wdfit, 'k') + f3, = plt.plot(checkedPpicks, checkedSPtimes, 'ko') + f4, = plt.plot(checkedPpicks, wdfit2, 'g') + plt.title('Wadati-Diagram, %d S-P Times, Vp/Vs(raw)=%5.2f,' \ + 'Vp/Vs(checked)=%5.2f' % (len(SPtimes), vpvsr, cvpvsr)) + plt.legend([f1, f2, f3, f4], ['Skipped S-Picks', 'Wadati 1', + 'Reliable S-Picks', 'Wadati 2'], loc='best') else: - plt.title('Wadati-Diagram, %d S-P Times' % len(SPtimes)) + plt.title('Wadati-Diagram, %d S-P Times' % len(SPtimes)) plt.ylabel('S-P Times [s]') plt.xlabel('P Times [s]') @@ -579,8 +582,8 @@ def wadaticheck(pickdic, dttolerance, iplot): def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot): ''' Function to detect spuriously picked noise peaks. - Uses RMS trace of all 3 components (if available) to determine, - how many samples [per cent] after P onset are below certain + Uses RMS trace of all 3 components (if available) to determine, + how many samples [per cent] after P onset are below certain threshold, calculated from noise level times noise factor. : param: X, time series (seismogram) @@ -612,7 +615,7 @@ def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot): print ("Checking signal length ...") if len(X) > 1: - # all three components available + # all three components available # make sure, all components have equal lengths ilen = min([len(X[0].data), len(X[1].data), len(X[2].data)]) x1 = X[0][0:ilen] @@ -639,7 +642,7 @@ def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot): numoverthr = len(np.where(rms[isignal] >= minsiglevel)[0]) if numoverthr >= minnum: - print ("checksignallength: Signal reached required length.") + print ("checksignallength: Signal reached required length.") returnflag = 1 else: print ("checksignallength: Signal shorter than required minimum signal length!") @@ -649,15 +652,15 @@ def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot): if iplot == 2: plt.figure(iplot) - p1, = plt.plot(t,rms, 'k') + p1, = plt.plot(t,rms, 'k') p2, = plt.plot(t[inoise], rms[inoise], 'c') p3, = plt.plot(t[isignal],rms[isignal], 'r') - p4, = plt.plot([t[isignal[0]], t[isignal[len(isignal)-1]]], \ - [minsiglevel, minsiglevel], 'g', linewidth=2) + p4, = plt.plot([t[isignal[0]], t[isignal[len(isignal)-1]]], + [minsiglevel, minsiglevel], 'g', linewidth=2) p5, = plt.plot([pick, pick], [min(rms), max(rms)], 'b', linewidth=2) - plt.legend([p1, p2, p3, p4, p5], ['RMS Data', 'RMS Noise Window', \ - 'RMS Signal Window', 'Minimum Signal Level', \ - 'Onset'], loc='best') + plt.legend([p1, p2, p3, p4, p5], ['RMS Data', 'RMS Noise Window', + 'RMS Signal Window', 'Minimum Signal Level', + 'Onset'], loc='best') plt.xlabel('Time [s] since %s' % X[0].stats.starttime) plt.ylabel('Counts') plt.title('Check for Signal Length, Station %s' % X[0].stats.station) @@ -729,32 +732,32 @@ def checkPonsets(pickdic, dttolerance, iplot): badjkmarker = 'badjkcheck' for i in range(0, len(goodstations)): # mark P onset as checked and keep P weight - pickdic[goodstations[i]]['P']['marked'] = goodmarker + pickdic[goodstations[i]]['P']['marked'] = goodmarker for i in range(0, len(badstations)): - # mark P onset and downgrade P weight to 9 - # (not used anymore) - pickdic[badstations[i]]['P']['marked'] = badmarker - pickdic[badstations[i]]['P']['weight'] = 9 + # mark P onset and downgrade P weight to 9 + # (not used anymore) + pickdic[badstations[i]]['P']['marked'] = badmarker + pickdic[badstations[i]]['P']['weight'] = 9 for i in range(0, len(badjkstations)): - # mark P onset and downgrade P weight to 9 - # (not used anymore) - pickdic[badjkstations[i]]['P']['marked'] = badjkmarker - pickdic[badjkstations[i]]['P']['weight'] = 9 + # mark P onset and downgrade P weight to 9 + # (not used anymore) + pickdic[badjkstations[i]]['P']['marked'] = badjkmarker + pickdic[badjkstations[i]]['P']['weight'] = 9 checkedonsets = pickdic if iplot > 1: - p1, = plt.plot(np.arange(0, len(Ppicks)), Ppicks, 'r+', markersize=14) + p1, = plt.plot(np.arange(0, len(Ppicks)), Ppicks, 'r+', markersize=14) p2, = plt.plot(igood, np.array(Ppicks)[igood], 'g*', markersize=14) - p3, = plt.plot([0, len(Ppicks) - 1], [pmedian, pmedian], 'g', \ - linewidth=2) + p3, = plt.plot([0, len(Ppicks) - 1], [pmedian, pmedian], 'g', + linewidth=2) for i in range(0, len(Ppicks)): - plt.text(i, Ppicks[i] + 0.2, stations[i]) + plt.text(i, Ppicks[i] + 0.2, stations[i]) plt.xlabel('Number of P Picks') plt.ylabel('Onset Time [s] from 1.1.1970') - plt.legend([p1, p2, p3], ['Skipped P Picks', 'Good P Picks', 'Median'], \ - loc='best') + plt.legend([p1, p2, p3], ['Skipped P Picks', 'Good P Picks', 'Median'], + loc='best') plt.title('Check P Onsets') plt.show() raw_input() @@ -789,37 +792,37 @@ def jackknife(X, phi, h): g = len(X) / h if type(g) is not int: - print ("jackknife: Cannot divide quantity X in equal sized subgroups!") + print ("jackknife: Cannot divide quantity X in equal sized subgroups!") print ("Choose another size for subgroups!") return PHI_jack, PHI_pseudo, PHI_sub else: - # estimator of undisturbed spot check - if phi == 'MEA': - phi_sc = np.mean(X) + # estimator of undisturbed spot check + if phi == 'MEA': + phi_sc = np.mean(X) elif phi == 'VAR': - phi_sc = np.var(X) + phi_sc = np.var(X) elif phi == 'MED': - phi_sc = np.median(X) + phi_sc = np.median(X) - # estimators of subgroups + # estimators of subgroups PHI_pseudo = [] PHI_sub = [] for i in range(0, g - 1): - # subgroup i, remove i-th sample - xx = X[:] - del xx[i] - # calculate estimators of disturbed spot check - if phi == 'MEA': - phi_sub = np.mean(xx) - elif phi == 'VAR': - phi_sub = np.var(xx) - elif phi == 'MED': - phi_sub = np.median(xx) + # subgroup i, remove i-th sample + xx = X[:] + del xx[i] + # calculate estimators of disturbed spot check + if phi == 'MEA': + phi_sub = np.mean(xx) + elif phi == 'VAR': + phi_sub = np.var(xx) + elif phi == 'MED': + phi_sub = np.median(xx) - PHI_sub.append(phi_sub) - # pseudo values - phi_pseudo = g * phi_sc - ((g - 1) * phi_sub) - PHI_pseudo.append(phi_pseudo) + PHI_sub.append(phi_sub) + # pseudo values + phi_pseudo = g * phi_sc - ((g - 1) * phi_sub) + PHI_pseudo.append(phi_pseudo) # jackknife estimator PHI_jack = np.mean(PHI_pseudo) @@ -899,29 +902,29 @@ def checkZ4S(X, pick, zfac, checkwin, iplot): # vertical P-coda level must exceed horizontal P-coda level # zfac times encodalevel if zcodalevel < minsiglevel: - print ("checkZ4S: Maybe S onset? Skip this P pick!") + print ("checkZ4S: Maybe S onset? Skip this P pick!") else: print ("checkZ4S: P onset passes checkZ4S test!") returnflag = 1 if iplot > 1: - te = np.arange(0, edat[0].stats.npts / edat[0].stats.sampling_rate, + te = np.arange(0, edat[0].stats.npts / edat[0].stats.sampling_rate, edat[0].stats.delta) - tn = np.arange(0, ndat[0].stats.npts / ndat[0].stats.sampling_rate, + tn = np.arange(0, ndat[0].stats.npts / ndat[0].stats.sampling_rate, ndat[0].stats.delta) - plt.plot(tz, z / max(z), 'k') + plt.plot(tz, z / max(z), 'k') plt.plot(tz[isignal], z[isignal] / max(z), 'r') plt.plot(te, edat[0].data / max(edat[0].data) + 1, 'k') plt.plot(te[isignal], edat[0].data[isignal] / max(edat[0].data) + 1, 'r') plt.plot(tn, ndat[0].data / max(ndat[0].data) + 2, 'k') plt.plot(tn[isignal], ndat[0].data[isignal] / max(ndat[0].data) + 2, 'r') - plt.plot([tz[isignal[0]], tz[isignal[len(isignal) - 1]]], \ - [minsiglevel / max(z), minsiglevel / max(z)], 'g', \ - linewidth=2) + plt.plot([tz[isignal[0]], tz[isignal[len(isignal) - 1]]], + [minsiglevel / max(z), minsiglevel / max(z)], 'g', + linewidth=2) plt.xlabel('Time [s] since %s' % zdat[0].stats.starttime) plt.ylabel('Normalized Counts') - plt.yticks([0, 1, 2], [zdat[0].stats.channel, edat[0].stats.channel, \ - ndat[0].stats.channel]) + plt.yticks([0, 1, 2], [zdat[0].stats.channel, edat[0].stats.channel, + ndat[0].stats.channel]) plt.title('CheckZ4S, Station %s' % zdat[0].stats.station) plt.show() raw_input() diff --git a/pylot/core/read/__init__.py b/pylot/core/read/__init__.py index 8b137891..40a96afc 100644 --- a/pylot/core/read/__init__.py +++ b/pylot/core/read/__init__.py @@ -1 +1 @@ - +# -*- coding: utf-8 -*- diff --git a/pylot/core/read/inputs.py b/pylot/core/read/inputs.py index 1c824ea1..01e7cd3b 100644 --- a/pylot/core/read/inputs.py +++ b/pylot/core/read/inputs.py @@ -208,8 +208,7 @@ class FilterOptions(object): def parseFilterOptions(self): if self.getFilterType(): - robject = {'type':self.getFilterType()} - robject['corners'] = self.getOrder() + robject = {'type': self.getFilterType(), 'corners': self.getOrder()} if len(self.getFreq()) > 1: robject['freqmin'] = self.getFreq()[0] robject['freqmax'] = self.getFreq()[1] diff --git a/pylot/core/read/io.py b/pylot/core/read/io.py index 926e3ca3..675aa150 100644 --- a/pylot/core/read/io.py +++ b/pylot/core/read/io.py @@ -73,8 +73,8 @@ def readPILOTEvent(phasfn=None, locfn=None, authority_id=None, **kwargs): stations = [stat for stat in phases['stat'][0:-1:3]] - event = createEvent(eventDate, loccinfo, None, 'earthquake', eventNum, - authority_id) + event = createEvent(eventDate, loccinfo, etype='earthquake', resID=eventNum, + authority_id=authority_id) lat = float(loc['LAT']) lon = float(loc['LON']) @@ -130,7 +130,7 @@ def readPILOTEvent(phasfn=None, locfn=None, authority_id=None, **kwargs): event.magnitudes.append(magnitude) return event - except AttributeError, e: + except AttributeError as e: raise AttributeError('{0} - Matlab LOC files {1} and {2} contains \ insufficient data!'.format(e, phasfn, locfn)) diff --git a/pylot/core/util/__init__.py b/pylot/core/util/__init__.py index 2d4a37e4..128de997 100755 --- a/pylot/core/util/__init__.py +++ b/pylot/core/util/__init__.py @@ -1 +1,2 @@ +# -*- coding: utf-8 -*- from pylot.core.util.version import get_git_version as _getVersionString diff --git a/pylot/core/util/testUtils.py b/pylot/core/util/testUtils.py index 9913c940..07d64925 100644 --- a/pylot/core/util/testUtils.py +++ b/pylot/core/util/testUtils.py @@ -1,3 +1,4 @@ +# -*- coding: utf-8 -*- ''' Created on 10.11.2014 @@ -23,4 +24,4 @@ class Test(unittest.TestCase): if __name__ == "__main__": #import sys;sys.argv = ['', 'Test.testName'] - unittest.main() \ No newline at end of file + unittest.main() diff --git a/pylot/core/util/testWidgets.py b/pylot/core/util/testWidgets.py index f5582362..b20190e1 100644 --- a/pylot/core/util/testWidgets.py +++ b/pylot/core/util/testWidgets.py @@ -1,3 +1,4 @@ +# -*- coding: utf-8 -*- ''' Created on 10.11.2014 @@ -15,4 +16,4 @@ class Test(unittest.TestCase): if __name__ == "__main__": #import sys;sys.argv = ['', 'Test.testName'] - unittest.main() \ No newline at end of file + unittest.main() diff --git a/pylot/core/util/thread.py b/pylot/core/util/thread.py index 6b0342f7..7ce1513e 100644 --- a/pylot/core/util/thread.py +++ b/pylot/core/util/thread.py @@ -1,3 +1,4 @@ +# -*- coding: utf-8 -*- import sys from PySide.QtCore import QThread, Signal diff --git a/pylot/core/util/utils.py b/pylot/core/util/utils.py index 2d3b0e8b..f5c9b5fd 100644 --- a/pylot/core/util/utils.py +++ b/pylot/core/util/utils.py @@ -1,4 +1,5 @@ #!/usr/bin/env python +# -*- coding: utf-8 -*- # # -*- coding: utf-8 -*- diff --git a/pylot/core/util/widgets.py b/pylot/core/util/widgets.py index 6e680d81..1a8ef94c 100644 --- a/pylot/core/util/widgets.py +++ b/pylot/core/util/widgets.py @@ -814,9 +814,10 @@ class PropertiesDlg(QDialog): if values is not None: self.setValues(values) - def setValues(self, tabValues): + @staticmethod + def setValues(tabValues): settings = QSettings() - for setting, value in tabValues.iteritems(): + for setting, value in tabValues.items(): settings.setValue(setting, value) settings.sync() diff --git a/testHelpForm.py b/testHelpForm.py index fe1d5a3e..2da58745 100755 --- a/testHelpForm.py +++ b/testHelpForm.py @@ -1,4 +1,5 @@ #!/usr/bin/env python +# -*- coding: utf-8 -*- import sys, time from PySide.QtGui import QApplication diff --git a/testPickDlg.py b/testPickDlg.py index e49f7675..5f8979b7 100755 --- a/testPickDlg.py +++ b/testPickDlg.py @@ -1,4 +1,5 @@ #!/usr/bin/env python +# -*- coding: utf-8 -*- import sys import matplotlib diff --git a/testPropDlg.py b/testPropDlg.py index d4c68871..b77de3b6 100755 --- a/testPropDlg.py +++ b/testPropDlg.py @@ -1,4 +1,5 @@ #!/usr/bin/env python +# -*- coding: utf-8 -*- import sys, time from PySide.QtGui import QApplication @@ -8,4 +9,4 @@ app = QApplication(sys.argv) win = PropertiesDlg() win.show() -app.exec_() \ No newline at end of file +app.exec_() diff --git a/testUIcomponents.py b/testUIcomponents.py index f422df81..920e9fd4 100755 --- a/testUIcomponents.py +++ b/testUIcomponents.py @@ -1,4 +1,6 @@ #!/usr/bin/env python +# -*- coding: utf-8 -*- + import sys, time from PySide.QtGui import QApplication @@ -9,7 +11,7 @@ dialogs = [FilterOptionsDialog, PropertiesDlg, HelpForm] app = QApplication(sys.argv) for dlg in dialogs: - win = dlg() - win.show() - time.sleep(1) - win.destroy() \ No newline at end of file + win = dlg() + win.show() + time.sleep(1) + win.destroy()