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

pylot/core/active/activeSeismoPick.py

@@ -135,26 +135,45 @@ 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])
+        
+        labelm = 'manual picks'
+        labela = 'automatic picks'
 
-        label = 'Difference to automatic picks [s]'
         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, 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])
+        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):
         '''
@@ -403,9 +422,9 @@ 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)
+                    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)

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

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,14 +358,18 @@ 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())
         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
 
@@ -373,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
 
@@ -663,7 +673,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 +698,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 +757,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 +785,10 @@ class SeismicShot(object):
                     [ax.get_ylim()[0],
                      ax.get_ylim()[1]],
                     'b:', label = 'latest')
-        ax.plot([0, self.getPick(traceID)],
+        if folm is not None:
-                 [folm * max(hoscf.getCF()), folm * max(hoscf.getCF())],
+            ax.plot([0, self.getPick(traceID)],
-                 'm:', label = 'folm = %s' %folm)
+                    [folm * max(hoscf.getCF()), folm * max(hoscf.getCF())],
+                    'm:', label = 'folm = %s' %folm)
         ax.set_xlabel('Time [s]')
         ax.legend()
 
@@ -808,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)
+        xaxis = np.arange(min(x) + step, max(x), step)
-        yaxis = np.arange(min(y), max(y), 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)
 
@@ -834,7 +847,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
 
@@ -886,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'
@@ -899,7 +912,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

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,23 +33,28 @@ 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, snrs, 'b.', markersize = 2.0, label = 'SNR values')
-    plt.plot(dists, snrthresholds, 'r.', markersize = 1, label = 'Fitted threshold')
+    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()
 
-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)
@@ -62,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):
     '''