Merge using remote to resolve conflicts

pylot/core/active/seismicshot.py

@@ -44,14 +44,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 +59,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:
@@ -147,10 +147,10 @@ class SeismicShot(object):
 
     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 +172,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]
 
@@ -257,16 +257,20 @@ class SeismicShot(object):
         else:
             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 ##########
+    def pickTraces(self, traceID, pickmethod, windowsize, folm = 0.6, HosAic = 'hos'): ########## input variables ##########
+        # LOCALMAX NOT IMPLEMENTED!
         '''
         Intitiate picking for a trace.
 
         :param: traceID
         :type: int
 
+        :param: pickmethod, use either 'threshold' or 'localmax' method. (localmax not yet implemented 04_08_15)
+        :type: string
+
         :param: cutwindow (equals HOScf 'cut' variable)
         :type: tuple
 
@@ -290,7 +294,13 @@ class SeismicShot(object):
 
         self.timeArray[traceID] = hoscf.getTimeArray()
 
-        aiccftime, hoscftime = self.threshold(hoscf, aiccf, windowsize, self.getPickwindow(traceID), folm)
+        if pickmethod == 'threshold':
+            aiccftime, hoscftime = self.threshold(hoscf, aiccf, windowsize, self.getPickwindow(traceID), folm)
+
+        #setpick = {'threshold':self.threshold,
+         #          'localmax':self.localmax}
+
+        #aiccftime, hoscftime = setpick[pickmethod](hoscf, aiccf, windowsize, pickwindow)
 
         setHosAic = {'hos': hoscftime,
                      'aic': aiccftime}
@@ -303,15 +313,15 @@ class SeismicShot(object):
         tsignal = self.getTsignal()
         tnoise = self.getPick(traceID) - tgap
 
-        (self.earliest[traceID], self.latest[traceID], tmp) = earllatepicker(self.getSingleStream(traceID), 
+        (self.earliest[traceID], self.latest[traceID], tmp) = earllatepicker(self.getSingleStream(traceID),
-                                                                             nfac, (tnoise, tgap, tsignal), 
+                                                                             nfac, (tnoise, tgap, tsignal),
                                                                              self.getPick(traceID))
 
     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 +347,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 +385,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 +399,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 +422,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
 
@@ -437,7 +447,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,15 +462,15 @@ 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
         # ++++++++++++++ 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()
@@ -474,12 +484,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 +520,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 +528,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
@@ -573,14 +584,14 @@ class SeismicShot(object):
         import matplotlib.pyplot as plt
         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')
-        
+
         hoscf = self.getHOScf(traceID)
         aiccf = self.getAICcf(traceID)
 
@@ -588,16 +599,16 @@ class SeismicShot(object):
         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)], 
+        ax1.plot([self.getPick(traceID), self.getPick(traceID)],
-                 [min(stream[0].data), 
+                 [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)], 
+        ax2.plot([self.getPick(traceID), self.getPick(traceID)],
-                 [min(np.minimum(hoscf.getCF(), aiccf.getCF())), 
+                 [min(np.minimum(hoscf.getCF(), aiccf.getCF())),
                   max(np.maximum(hoscf.getCF(), aiccf.getCF()))],
                  'r', label = 'mostlikely')
         ax2.plot([0, self.getPick(traceID)],
@@ -605,7 +616,7 @@ class SeismicShot(object):
                  'm:', label = 'folm = %s' %folm)
         plt.xlabel('Time [s]')
         plt.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.
@@ -644,7 +655,7 @@ class SeismicShot(object):
         if ax == None:
             fig = plt.figure()
             ax = plt.axes(projection = '3d')
-        
+
         xsrc, ysrc, zsrc = self.getSrcLoc()
 
         if contour == True:
@@ -656,12 +667,12 @@ 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):
         '''
         Plots a 2D matrix of the interpolated traveltimes. This needs less performance than plot3dttc
@@ -701,7 +712,7 @@ class SeismicShot(object):
             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:

pylot/core/active/surveyUtils.py

@@ -1,19 +1,69 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 
-import numpy as np
+from pylab import *
+startpos = []
+endpos = []
+
+def generateSurvey(obsdir, shotlist):
+    from obspy.core import read
+    from pylot.core.active import seismicshot
+
+    shot_dict = {}
+    for shotnumber in shotlist:       # loop over data files
+        # generate filenames and read manual picks to a list
+        obsfile = obsdir + str(shotnumber) + '_pickle.dat'
+        #obsfile = obsdir + str(shotnumber) + '.dat'
+
+        if not obsfile in shot_dict.keys():
+            shot_dict[shotnumber] = []
+        shot_dict[shotnumber] = seismicshot.SeismicShot(obsfile)
+        shot_dict[shotnumber].setParameters('shotnumber', shotnumber)
+
+    return shot_dict
+
+def setParametersForShots(cutwindow, tmovwind, tsignal, tgap, receiverfile, sourcefile, shot_dict):
+    for shot in shot_dict.values():
+        shot.setCut(cutwindow)
+        shot.setTmovwind(tmovwind)
+        shot.setTsignal(tsignal)
+        shot.setTgap(tgap)
+        shot.setRecfile(receiverfile)
+        shot.setSourcefile(sourcefile)
+        shot.setOrder(order = 4)
+
+def removeEmptyTraces(shot_dict):
+    filename = 'removeEmptyTraces.out'
+    filename2 = 'updateTraces.out'
+    outfile = open(filename, 'w')
+    outfile2 = open(filename2, 'w')
+    for shot in shot_dict.values():
+        del_traceIDs = shot.updateTraceList()
+        removed = shot.removeEmptyTraces()
+        if removed is not None:
+            outfile.writelines('shot: %s, removed empty traces: %s\n' %(shot.getShotnumber(), removed))
+            outfile2.writelines('shot: %s, removed traceID(s) %s because they were not found in the corresponding stream\n' %(shot.getShotnumber(), del_traceIDs))
+    print '\nremoveEmptyTraces, updateTraces: Finished! See %s and %s for more information of removed traces.\n' %(filename, filename2)
+    outfile.close()
+    outfile2.close()
+
 
 def readParameters(parfile, parameter):
     from ConfigParser import ConfigParser
     parameterConfig = ConfigParser()
     parameterConfig.read('parfile')
-    
+
-    value = parameterConfig.get('vars', parameter).split('#')[0]
+    value = parameterConfig.get('vars', parameter).split('\t')[0]
-    value = value.replace(" ", "")
 
     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 = []
@@ -34,6 +84,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 +96,98 @@ 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 linearInterp(dist_med, dist_start
+
+def exportFMTOMO(shot_dict, directory = 'FMTOMO_export', sourcefile = 'input_sf.in', ttFileExtension = '.tt'):
+    count = 0
+    fmtomo_factor = 1000 # transforming [m/s] -> [km/s]
+    LatAll = []; LonAll = []; DepthAll = []
+    srcfile = open(directory + '/' + sourcefile, 'w')
+    srcfile.writelines('%10s\n' %len(shot_dict)) # number of sources
+    for shotnumber in getShotlist(shot_dict):
+        shot = getShotForShotnumber(shot_dict, shotnumber)
+        ttfilename = str(shotnumber) + ttFileExtension
+        (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 %10s %10s\n' %(1, 1, ttfilename))
+        ttfile = open(directory + '/' + ttfilename, 'w')
+        traceIDlist = shot.getTraceIDlist()
+        traceIDlist.sort()
+        ttfile.writelines(str(countPickedTraces(shot)) + '\n')
+        for traceID in traceIDlist:
+            if shot.getPick(traceID) is not None:
+                pick = shot.getPick(traceID) * fmtomo_factor
+                delta = shot.getPickError(traceID) * fmtomo_factor
+                (x, y, z) = shot.getRecLoc(traceID)
+                ttfile.writelines('%20s %20s %20s %10s %10s\n' %(getAngle(y), getAngle(x), (-1)*z, pick, delta))
+                LatAll.append(getAngle(y)); LonAll.append(getAngle(x)); DepthAll.append((-1)*z)
+                count += 1
+        ttfile.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)'%(
+        min(DepthAll), max(DepthAll), min(LatAll), max(LatAll), min(LonAll), max(LonAll))
+
+def getShotlist(shot_dict):
+    shotlist = []
+    for shot in shot_dict.values():
+        shotlist.append(shot.getShotnumber())
+    shotlist.sort()
+    return shotlist
+
+def getShotForShotnumber(shot_dict, shotnumber):
+    for shot in shot_dict.values():
+        if shot.getShotnumber() == shotnumber:
+            return shot
+
+def getAngle(distance):
+    '''
+    Function returns the angle on a Sphere of the radius R = 6371 [km] for a distance [km].
+    '''
+    import numpy as np
+    PI = np.pi
+    R = 6371.
+    angle = distance * 180 / (PI * R)
+    return angle
+
+def countPickedTraces(shot):
+    numtraces = 0
+    for traceID in shot.getTraceIDlist():
+        if shot.getPick(traceID) is not None:
+            numtraces += 1
+    print "countPickedTraces: Found %s picked traces in shot number %s"  %(numtraces, shot.getShotnumber())
+    return numtraces
+
+def countAllPickedTraces(shot_dict):
+    traces = 0
+    for shot in shot_dict.values():
+        traces += countPickedTraces(shot)
+    return traces
 
 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():
@@ -75,3 +199,6 @@ def findTracesInRanges(shot_dict, distancebin, pickbin):
                     shots_found[shot.getShotnumber()].append(traceID)
 
     return shots_found
+
+
+