From 0b1f16866bb45859bebb651343a1f448009bd621 Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Fri, 11 Dec 2015 11:08:06 +0100 Subject: [PATCH 1/5] added plotHist --- pylot/core/active/activeSeismoPick.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/pylot/core/active/activeSeismoPick.py b/pylot/core/active/activeSeismoPick.py index 3cf59e80..ce8ea46d 100644 --- a/pylot/core/active/activeSeismoPick.py +++ b/pylot/core/active/activeSeismoPick.py @@ -155,6 +155,17 @@ class Survey(object): ax.set_ylabel('Time [s]') ax.text(0.5, 0.95, 'Plot of all MANUAL picks', transform=ax.transAxes, horizontalalignment='center') + def plotHist(self): + import matplotlib.pyplot as plt + plt.interactive(True) + diffs = [] + for shot in self.data.values(): + for traceID in shot.getTraceIDlist(): + if shot.getPickFlag(traceID) == 1 and shot.getManualPickFlag(traceID) == 1: + diffs.append(self.getDiffsFromManual()[shot][traceID]) + plt.hist(diffs, 20) + plt.title('Histogram of the differences between automatic and manual pick') + plt.xlabel('Difference in time (auto - manual) [s]') def pickAllShots(self, windowsize, HosAic = 'hos', vmin = 333, vmax = 5500, folm = 0.6): ''' From 3b4e1dcd1ea4c6ef2eb97b3d1428d7eb9516e8ce Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Mon, 4 Jan 2016 15:46:24 +0100 Subject: [PATCH 2/5] added changes for manual picks --- pylot/core/active/activeSeismoPick.py | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/pylot/core/active/activeSeismoPick.py b/pylot/core/active/activeSeismoPick.py index ce8ea46d..5ae1cfa9 100644 --- a/pylot/core/active/activeSeismoPick.py +++ b/pylot/core/active/activeSeismoPick.py @@ -135,35 +135,41 @@ class Survey(object): def plotDiffs(self): import matplotlib.pyplot as plt - diffs = []; dists = []; picks = [] + diffs = []; dists = []; mpicks = []; picks = [] diffsDic = self.getDiffsFromManual() for shot in self.data.values(): for traceID in shot.getTraceIDlist(): if shot.getPickFlag(traceID) == 1 and shot.getManualPickFlag(traceID) == 1: dists.append(shot.getDistance(traceID)) - picks.append(shot.getManualPick(traceID)) + mpicks.append(shot.getManualPick(traceID)) + picks.append(shot.getPick(traceID)) diffs.append(diffsDic[shot][traceID]) - label = 'Difference to automatic picks [s]' + labelm = 'manual picks' + labela = 'automatic picks' fig = plt.figure() ax = fig.add_subplot(111) - sc = ax.scatter(dists, picks, c = diffs, s=5, edgecolors='none', label = label) + sc_a = ax.scatter(dists, picks, c = '0.5', s=10, edgecolors='none', label = labela, alpha = 0.3) + sc = ax.scatter(dists, mpicks, c = diffs, s=5, edgecolors='none', label = labelm) cbar = plt.colorbar(sc, fraction=0.05) - cbar.set_label(label) + cbar.set_label(labelm) ax.set_xlabel('Distance [m]') ax.set_ylabel('Time [s]') ax.text(0.5, 0.95, 'Plot of all MANUAL picks', transform=ax.transAxes, horizontalalignment='center') - def plotHist(self): + def plotHist(self, nbins = 20, ax = None): import matplotlib.pyplot as plt plt.interactive(True) diffs = [] + if ax == None: + fig = plt.figure() + ax = fig.add_subplot(111) for shot in self.data.values(): for traceID in shot.getTraceIDlist(): if shot.getPickFlag(traceID) == 1 and shot.getManualPickFlag(traceID) == 1: diffs.append(self.getDiffsFromManual()[shot][traceID]) - plt.hist(diffs, 20) + plt.hist(diffs, nbins) plt.title('Histogram of the differences between automatic and manual pick') plt.xlabel('Difference in time (auto - manual) [s]') @@ -407,7 +413,7 @@ class Survey(object): #shot_dict[shotnumber].plot3dttc(ax = ax, plotpicks = True) ax = fig.add_subplot(3, 4, index) if mode == '3d': - self.getShot(shotnumber).matshow(ax = ax, colorbar = False, annotations = True) + self.getShot(shotnumber).matshow(ax = ax, colorbar = False, annotations = True, legend = False) elif mode == '2d': self.getShot(shotnumber).plot2dttc(ax) self.getShot(shotnumber).plotmanual2dttc(ax) From 3c1be950b94b63858f5e80dd991ad00d8d1ac0a8 Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Mon, 4 Jan 2016 15:47:05 +0100 Subject: [PATCH 3/5] removed most of the folm = 0.6 default values --- pylot/core/active/seismicshot.py | 29 ++++++++++++++++++----------- 1 file changed, 18 insertions(+), 11 deletions(-) diff --git a/pylot/core/active/seismicshot.py b/pylot/core/active/seismicshot.py index 2f6595f9..0bed77d4 100644 --- a/pylot/core/active/seismicshot.py +++ b/pylot/core/active/seismicshot.py @@ -37,6 +37,7 @@ class SeismicShot(object): self.traces4plot = {} self.paras = {} self.paras['shotname'] = obsfile + self.folm = None def removeEmptyTraces(self): traceIDs = [] @@ -279,7 +280,7 @@ class SeismicShot(object): #raise ValueError('ambigious or empty traceID: %s' % traceID) - def pickTraces(self, traceID, windowsize, folm = 0.6, HosAic = 'hos'): ########## input variables ########## + def pickTraces(self, traceID, windowsize, folm, HosAic = 'hos'): ########## input variables ########## # LOCALMAX NOT IMPLEMENTED! ''' Intitiate picking for a trace. @@ -299,7 +300,7 @@ class SeismicShot(object): :param: windowsize, window around the returned HOS picktime, to search for the AIC minumum :type: 'tuple' - :param: folm, fraction of local maximumm (default = 0.6) + :param: folm, fraction of local maximumm :type: 'real' :param: HosAic, get hos or aic pick (can be 'hos'(default) or 'aic') @@ -308,6 +309,8 @@ class SeismicShot(object): hoscf = self.getHOScf(traceID) ### determination of both, HOS and AIC (need to change threshold-picker) ### aiccf = self.getAICcf(traceID) + self.folm = folm + self.timeArray[traceID] = hoscf.getTimeArray() aiccftime, hoscftime = self.threshold(hoscf, aiccf, windowsize, self.getPickwindow(traceID), folm) setHosAic = {'hos': hoscftime, @@ -335,7 +338,7 @@ class SeismicShot(object): # self.picks[traceID]['spe'] *= 0.5 # TEST OF 1/2 PICKERROR - def threshold(self, hoscf, aiccf, windowsize, pickwindow, folm = 0.6): + def threshold(self, hoscf, aiccf, windowsize, pickwindow, folm): ''' 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. @@ -355,7 +358,7 @@ class SeismicShot(object): :param: cutwindow [seconds], cut a part of the trace as in Characteristic Function :type: 'tuple' - :param: folm, fraction of local maximum (default = 0.6) + :param: folm, fraction of local maximum :type: 'real' ''' hoscflist = list(hoscf.getCF()) @@ -663,7 +666,7 @@ class SeismicShot(object): ax.legend() ax.text(0.05, 0.9, 'SNR: %s' %snr, transform = ax.transAxes) - def plot_traces(self, traceID, folm = 0.6): ########## 2D, muss noch mehr verbessert werden ########## + def plot_traces(self, traceID): ########## 2D, muss noch mehr verbessert werden ########## from matplotlib.widgets import Button def onclick(event): @@ -688,6 +691,8 @@ class SeismicShot(object): def cleanup(event): self.traces4plot[traceID] = {} + folm = self.folm + fig = plt.figure() ax1 = fig.add_subplot(2,1,1) ax2 = fig.add_subplot(2,1,2, sharex = ax1) @@ -745,7 +750,7 @@ class SeismicShot(object): ax.legend() return ax - def _drawCFs(self, traceID, folm, refresh = False): + def _drawCFs(self, traceID, folm = None, refresh = False): hoscf = self.getHOScf(traceID) aiccf = self.getAICcf(traceID) ax = self.traces4plot[traceID]['ax2'] @@ -773,9 +778,10 @@ class SeismicShot(object): [ax.get_ylim()[0], ax.get_ylim()[1]], 'b:', label = 'latest') - ax.plot([0, self.getPick(traceID)], - [folm * max(hoscf.getCF()), folm * max(hoscf.getCF())], - 'm:', label = 'folm = %s' %folm) + if folm is not None: + ax.plot([0, self.getPick(traceID)], + [folm * max(hoscf.getCF()), folm * max(hoscf.getCF())], + 'm:', label = 'folm = %s' %folm) ax.set_xlabel('Time [s]') ax.legend() @@ -834,7 +840,7 @@ class SeismicShot(object): plotmethod[method](*args) - def matshow(self, ax = None, step = 0.5, method = 'linear', plotRec = True, annotations = True, colorbar = True): + def matshow(self, ax = None, step = 0.5, method = 'linear', plotRec = True, annotations = True, colorbar = True, legend = True): ''' Plots a 2D matrix of the interpolated traveltimes. This needs less performance than plot3dttc @@ -899,7 +905,8 @@ class SeismicShot(object): cbar = plt.colorbar(sc) cbar.set_label('Time [s]') - ax.legend() + if legend == True: + ax.legend() ax.set_xlabel('X') ax.set_ylabel('Y') ax.plot(self.getSrcLoc()[0], self.getSrcLoc()[1],'*k', markersize = 15) # plot source location From acd9e942f97cbd3618fab9f5c03bbbc37f2acd3a Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Mon, 4 Jan 2016 15:47:46 +0100 Subject: [PATCH 4/5] improved SNR plots --- pylot/core/active/surveyUtils.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/pylot/core/active/surveyUtils.py b/pylot/core/active/surveyUtils.py index 12f5cc39..a3081517 100644 --- a/pylot/core/active/surveyUtils.py +++ b/pylot/core/active/surveyUtils.py @@ -16,11 +16,13 @@ def setArtificialPick(shot_dict, traceID, pick): def fitSNR4dist(shot_dict, shiftdist = 30, shiftSNR = 100): import numpy as np + import matplotlib.pyplot as plt dists = [] picks = [] snrs = [] snr_sqrt_inv = [] snrthresholds = [] + snrBestFit = [] for shot in shot_dict.values(): for traceID in shot.getTraceIDlist(): if shot.getSNR(traceID)[0] >= 1: @@ -31,18 +33,23 @@ def fitSNR4dist(shot_dict, shiftdist = 30, shiftSNR = 100): fit = np.polyfit(dists, snr_sqrt_inv, 1) fit_fn = np.poly1d(fit) for dist in dists: + snrBestFit.append((1/(fit_fn(dist)**2))) dist += shiftdist snrthresholds.append((1/(fit_fn(dist)**2)) - shiftSNR * np.exp(-0.05 * dist)) - plotFittedSNR(dists, snrthresholds, snrs) + plotFittedSNR(dists, snrthresholds, snrs, snrBestFit) return fit_fn #### ZU VERBESSERN, sollte fertige funktion wiedergeben -def plotFittedSNR(dists, snrthresholds, snrs): +def plotFittedSNR(dists, snrthresholds, snrs, snrBestFit): import matplotlib.pyplot as plt plt.interactive(True) fig = plt.figure() - plt.plot(dists, snrs, '.', markersize = 1.0, label = 'SNR values') - plt.plot(dists, snrthresholds, 'r.', markersize = 1, label = 'Fitted threshold') + plt.plot(dists, snrs, 'b.', markersize = 2.0, label = 'SNR values') + dists.sort() + snrthresholds.sort(reverse = True) + snrBestFit.sort(reverse = True) + plt.plot(dists, snrthresholds, 'r', markersize = 1, label = 'Fitted threshold') + plt.plot(dists, snrBestFit, 'k', markersize = 1, label = 'Best fitted curve') plt.xlabel('Distance[m]') plt.ylabel('SNR') plt.legend() From ce705888f198e41c4fb8e6501a8146a3c3b55f6c Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Tue, 19 Jan 2016 10:31:06 +0100 Subject: [PATCH 5/5] bugfixes and other not further specified changes --- pylot/core/active/activeSeismoPick.py | 5 +++-- pylot/core/active/seismicArrayPreparation.py | 13 +++++++++---- pylot/core/active/seismicshot.py | 17 ++++++++++++----- pylot/core/active/surveyUtils.py | 10 ++++++++-- 4 files changed, 32 insertions(+), 13 deletions(-) diff --git a/pylot/core/active/activeSeismoPick.py b/pylot/core/active/activeSeismoPick.py index 5632664f..2e8c4058 100644 --- a/pylot/core/active/activeSeismoPick.py +++ b/pylot/core/active/activeSeismoPick.py @@ -170,9 +170,10 @@ class Survey(object): for traceID in shot.getTraceIDlist(): if shot.getPickFlag(traceID) == 1 and shot.getManualPickFlag(traceID) == 1: diffs.append(self.getDiffsFromManual()[shot][traceID]) - plt.hist(diffs, nbins) + hist = plt.hist(diffs, nbins, histtype = 'step', normed = True, stacked = True) plt.title('Histogram of the differences between automatic and manual pick') plt.xlabel('Difference in time (auto - manual) [s]') + return diffs def pickAllShots(self, windowsize, HosAic = 'hos', vmin = 333, vmax = 5500, folm = 0.6): ''' @@ -421,7 +422,7 @@ class Survey(object): #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) + ax = fig.add_subplot(rows, columns, index) if mode == '3d': self.getShot(shotnumber).matshow(ax = ax, colorbar = False, annotations = True, legend = False) elif mode == '2d': diff --git a/pylot/core/active/seismicArrayPreparation.py b/pylot/core/active/seismicArrayPreparation.py index 607d904f..22cd8486 100644 --- a/pylot/core/active/seismicArrayPreparation.py +++ b/pylot/core/active/seismicArrayPreparation.py @@ -747,6 +747,9 @@ class SeisArray(object): return surface def exportAll(self, filename = 'interpolated_receivers.out'): + ''' + Exports all receivers to an input file for ActiveSeismoPick3D. + ''' recfile_out = open(filename, 'w') count = 0 for traceID in self.getReceiverCoordinates().keys(): @@ -803,7 +806,7 @@ class SeisArray(object): xrc, yrc, zrc = self.getReceiverLists() xsc, ysc, zsc = self.getSourceLocsLists() - plt.title('3D plot of seismic array %s'%self.recfile) + plt.title('3D plot of seismic array.') if len(xmt) > 0: ax.plot(xmt, ymt, zmt, 'b.', markersize = 10, label = 'measured topo points') if len(xrc) > 0: @@ -812,7 +815,7 @@ class SeisArray(object): ax.plot(xmr, ymr, zmr, 'ro', label = 'measured receivers') if len(xsc) > 0: ax.plot(xsc, ysc, zsc, 'b*', label = 'shot locations') - ax.set_xlabel('X'); ax.set_ylabel('Y'); ax.set_zlabel('elevation') + ax.set_xlabel('X [m]'); ax.set_ylabel('Y [m]'); ax.set_zlabel('Z [m]') ax.legend() return ax @@ -842,13 +845,15 @@ class SeisArray(object): 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)) + surf = ax.plot_surface(xgrid, ygrid, zgrid, linewidth = 0, cmap = cm.jet, vmin = min(z), vmax = max(z)) + cbar = plt.colorbar(surf) + cbar.set_label('Elevation [m]') if exag == False: 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.set_xlabel('X [m]'); ax.set_ylabel('Y [m]'); ax.set_zlabel('Z [m]') ax.legend() return ax diff --git a/pylot/core/active/seismicshot.py b/pylot/core/active/seismicshot.py index 0bed77d4..8fbdbfb8 100644 --- a/pylot/core/active/seismicshot.py +++ b/pylot/core/active/seismicshot.py @@ -365,7 +365,11 @@ class SeismicShot(object): leftb = int(pickwindow[0] / self.getCut()[1] * len(hoscflist)) rightb = int(pickwindow[1] / self.getCut()[1] * len(hoscflist)) - threshold = folm * max(hoscflist[leftb : rightb]) # combination of local maximum and threshold + #threshold = folm * max(hoscflist[leftb : rightb]) # combination of local maximum and threshold + + ### TEST TEST + threshold = folm * (max(hoscflist[leftb : rightb]) - min(hoscflist[leftb : rightb])) + min(hoscflist[leftb : rightb]) # combination of local maximum and threshold + ### TEST TEST m = leftb @@ -376,7 +380,10 @@ class SeismicShot(object): lb = max(0, m - windowsize[0]) # if window exceeds t = 0 aiccfcut = list(aiccf.getCF())[lb : m + windowsize[1]] - n = aiccfcut.index(min(aiccfcut)) + if len(aiccfcut) > 0: + n = aiccfcut.index(min(aiccfcut)) + else: + n = 0 m = lb + n @@ -814,8 +821,8 @@ class SeismicShot(object): y.append(self.getRecLoc(traceID)[1]) z.append(self.getPick(traceID)) - xaxis = np.arange(min(x), max(x), step) - yaxis = np.arange(min(y), max(y), step) + xaxis = np.arange(min(x) + step, max(x), step) + yaxis = np.arange(min(y) + step, max(y), step) xgrid, ygrid = np.meshgrid(xaxis, yaxis) zgrid = griddata((x, y), z, (xgrid, ygrid), method = method) @@ -892,9 +899,9 @@ class SeismicShot(object): ax.text(0.5, 0.95, 'shot: %s' %self.getShotnumber(), transform = ax.transAxes , horizontalalignment = 'center') sc = ax.scatter(x, y, c = z, s = 30, label = 'picked shots', vmin = tmin, vmax = tmax, cmap = cmap, linewidths = 1.5) + label = None for xyz in zip(xcut, ycut, zcut): x, y, z = xyz - label = None if z > tmax: count += 1 z = 'w' diff --git a/pylot/core/active/surveyUtils.py b/pylot/core/active/surveyUtils.py index a3081517..c4180e10 100644 --- a/pylot/core/active/surveyUtils.py +++ b/pylot/core/active/surveyUtils.py @@ -54,7 +54,7 @@ def plotFittedSNR(dists, snrthresholds, snrs, snrBestFit): plt.ylabel('SNR') plt.legend() -def setFittedSNR(shot_dict, shiftdist = 30, shiftSNR = 100, p1 = 0.004, p2 = -0.0007): +def setDynamicFittedSNR(shot_dict, shiftdist = 30, shiftSNR = 100, p1 = 0.004, p2 = -0.0007): import numpy as np minSNR = 2.5 #fit_fn = fitSNR4dist(shot_dict) @@ -69,8 +69,14 @@ def setFittedSNR(shot_dict, shiftdist = 30, shiftSNR = 100, p1 = 0.004, p2 = -0. shot.setSNRthreshold(traceID, minSNR) else: shot.setSNRthreshold(traceID, snrthreshold) - print "setFittedSNR: Finished setting of fitted SNR-threshold" + print "setDynamicFittedSNR: Finished setting of fitted SNR-threshold" +def setConstantSNR(shot_dict, snrthreshold = 2.5): + import numpy as np + for shot in shot_dict.values(): + for traceID in shot.getTraceIDlist(): + shot.setSNRthreshold(traceID, snrthreshold) + print "setConstantSNR: Finished setting of SNR threshold to a constant value of %s"%snrthreshold def findTracesInRanges(shot_dict, distancebin, pickbin): '''