added new files to pylot active

pylot/core/active/activeSeismoPick.py

@@ -0,0 +1,218 @@
+from pylot.core.active import seismicshot
+import numpy as np
+
+class Survey(object):
+    def __init__(self, path, sourcefile, receiverfile, useDefaultParas = False):
+        '''
+        The Survey Class contains all shots [type: seismicshot] of a survey
+        as well as the aquisition geometry and the topography.
+        '''
+        self.data = {}
+        self._topography = None
+        self._recfile = receiverfile
+        self._sourcefile = sourcefile
+        self._obsdir = path
+        self._generateSurvey()
+        if useDefaultParas == True: 
+            self.setParametersForShots()
+        self._removeAllEmptyTraces()
+        self._updateShots()
+        
+    def _generateSurvey(self):
+        from obspy.core import read
+
+        shot_dict = {}
+        shotlist = self.getShotlist()
+        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():
+                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 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.")
+        # CHANGE this later. Parameters only needed for survey, not for each shot.
+        for shot in self.data.values():
+            shot.setCut(cutwindow)
+            shot.setTmovwind(tmovwind)
+            shot.setTsignal(tsignal)
+            shot.setTgap(tgap)
+            shot.setRecfile(self.getPath() + self.getReceiverfile())
+            shot.setSourcefile(self.getPath() + self.getSourcefile())
+            shot.setOrder(order = 4)
+        print ("setParametersForShots: Parameters set to:\n"
+               "cutwindow = %s, tMovingWindow = %f, tsignal = %f, tgap = %f"
+               %(cutwindow, tmovwind, tsignal, tgap))
+
+    def _removeAllEmptyTraces(self):
+        filename = 'removeEmptyTraces.out'
+        count = 0
+        for shot in self.data.values():
+            removed = shot.removeEmptyTraces()
+            if removed is not None:
+                if count == 0: outfile = open(filename, 'w')
+                count += 1
+                outfile.writelines('shot: %s, removed empty traces: %s\n' 
+                                   %(shot.getShotnumber(), removed))
+        print ("\nremoveEmptyTraces: Finished! Removed %d traces" %count)
+        if count > 0:
+            print ("See %s for more information "
+                   "on removed traces."%(filename))
+            outfile.close()
+
+    def _updateShots(self):
+        filename = 'updateShots.out'
+        count = 0; countTraces = 0
+        for shot in self.data.values():
+            del_traceIDs = shot.updateTraceList()
+            if len(del_traceIDs) > 0:
+                if count == 0: outfile = open(filename, 'w')
+                count += 1
+                countTraces += len(del_traceIDs)
+                outfile.writelines("shot: %s, removed traceID(s) %s because "
+                                   "they were not found in the corresponding stream\n" 
+                                   %(shot.getShotnumber(), del_traceIDs))
+
+        print ("\nupdateShots: Finished! Updated %d shots and removed "
+               "%d traces" %(count, countTraces))
+        if count > 0:
+            print ("See %s for more information "
+                   "on removed traces."%(filename))
+            outfile.close()
+
+    def pickAllShots(self, HosAic = 'hos', vmin = 333, vmax = 5500):
+        '''
+        Automatically pick all traces of all shots of the survey.
+        '''
+        from datetime import datetime
+        count = 0
+
+        for shot in self.data.values():
+            tstartpick = datetime.now(); count += 1
+            for traceID in shot.getTraceIDlist():
+                distance = shot.getDistance(traceID) # receive distance
+
+                pickwin_used = shot.getCut()
+                cutwindow = shot.getCut()
+
+                # 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, pickmethod, 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:
+                    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)
+
+
+
+
+    def setArtificialPick(self, traceID, pick):
+        for shot in self.data.values():
+            shot.setPick(traceID, pick)
+            shot.setPickwindow(traceID, shot.getCut())
+
+    def countAllTraces(self):
+        numtraces = 0
+        for line in self.getShotlist():
+            for line in self.getReceiverlist():
+                numtraces += 1
+        return numtraces
+
+    def getShotlist(self):
+        filename = self.getPath() + self.getSourcefile()
+        srcfile = open(filename,'r')
+        shotlist = []
+        for line in srcfile.readlines():
+            line = line.split()
+            shotlist.append(int(line[0]))
+
+        return shotlist
+
+    def getReceiverlist(self):
+        filename = self.getPath() + self.getReceiverfile()
+        recfile = open(filename,'r')
+        reclist = []
+        for line in recfile.readlines():
+            line = line.split()
+            reclist.append(int(line[0]))
+
+        return reclist
+
+    def getShotDict(self):
+        return self.data
+
+    def getShot(self, shotnumber):
+        return self.data[shotnumber]
+
+    def getSourcefile(self):
+        return self._sourcefile
+
+    def getReceiverfile(self):
+        return self._recfile
+    
+    def getPath(self):
+        return self._obsdir
+
+    def getStats(self):
+        info_dict = {}
+        for shot in self.data.values():
+            pickedTraces = 0
+            snrlist = []
+            dist = []
+            numtraces = len(shot.getTraceIDlist())
+            for traceID in shot.getTraceIDlist():
+                snrlist.append(shot.getSNR(traceID)[0])
+                dist.append(shot.getDistance(traceID))
+                if shot.getPick(traceID) is not None:
+                    pickedTraces += 1
+            info_dict[shot.getShotnumber()] = {'numtraces': numtraces,
+                                               'picked traces': [pickedTraces,
+                                                                 '%2.2f %%'%(pickedTraces /
+                                                                             numtraces * 100)],
+                                               'mean SNR': np.mean(snrlist),
+                                               'mean distance': np.mean(dist)}
+
+        return info_dict
+
+    def saveSurvey(self, filename = 'survey.pickle'):
+        import cPickle
+        outfile = open(filename, 'wb')
+
+        cPickle.dump(self, outfile, -1)
+        print('saved Survey to file %s'%(filename))
+        
+    @staticmethod    
+    def from_pickle(filename):
+        import cPickle
+        infile = open(filename, 'rb')
+        return cPickle.load(infile)

pylot/core/active/fmtomo2vtk.py

@@ -0,0 +1,204 @@
+def vgrids2VTK(inputfile = 'vgrids.in', outputfile = 'vgrids.vtk'):
+    '''
+    Generate a vtk-file readable by e.g. paraview from FMTOMO output vgrids.in
+    '''
+    def getDistance(angle):
+        PI = np.pi
+        R = 6371.
+        distance = angle / 180 * (PI * R)
+        return distance
+
+    def readNumberOfPoints(filename):
+        fin = open(filename, 'r')
+        vglines = fin.readlines()
+
+        nR = int(vglines[1].split()[0])
+        nTheta = int(vglines[1].split()[1])
+        nPhi = int(vglines[1].split()[2])
+
+        fin.close()
+        return nR, nTheta, nPhi
+
+    def readDelta(filename):
+        fin = open(filename, 'r')
+        vglines = fin.readlines()
+
+        dR = float(vglines[2].split()[0])
+        dTheta = float(vglines[2].split()[1])
+        dPhi = float(vglines[2].split()[2])
+
+        fin.close()
+        return dR, dTheta, dPhi
+
+    def readStartpoints(filename):
+        fin = open(filename, 'r')
+        vglines = fin.readlines()
+
+        sR = float(vglines[3].split()[0])
+        sTheta = float(vglines[3].split()[1])
+        sPhi = float(vglines[3].split()[2])
+
+        fin.close()
+        return sR, sTheta, sPhi
+
+    def readVelocity(filename):
+        vel = []; count = 0
+        fin = open(filename, 'r')
+        vglines = fin.readlines()
+
+        for line in vglines:
+            count += 1
+            if count > 4:
+                vel.append(float(line.split()[0]))
+
+        print("Read %d points out of file: %s" %(count - 4, filename))
+        return vel
+
+    R = 6371 # earth radius
+    outfile = open(outputfile, 'w')
+
+    # Theta, Phi in radians, R in km
+    nR, nTheta, nPhi = readNumberOfPoints(inputfile)
+    dR, dTheta, dPhi = readDelta(inputfile)
+    sR, sTheta, sPhi = readStartpoints(inputfile)
+    vel = readVelocity(inputfile)
+
+    nX = nPhi; nY = nTheta; nZ = nR
+
+    sZ = sR - R
+    sX = getDistance(np.rad2deg(sPhi))
+    sY = getDistance(np.rad2deg(sTheta))
+
+    dX = getDistance(np.rad2deg(dPhi))
+    dY = getDistance(np.rad2deg(dTheta))
+    
+    xGrid = np.linspace(sX, sX + (dX * nX), nX)
+    yGrid = np.linspace(sZ, sZ + (nY * dY), nY)
+    zGrid = np.linspace(sZ, sZ + (nR * dR), nR)
+    nPoints = len(xGrid) * len(yGrid) * len(zGrid)
+
+    # write header
+    print("Writing header for VTK file...")
+    outfile.writelines('# vtk DataFile Version 3.1\n')
+    outfile.writelines('Velocity on FMTOMO vgrids.in points\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 z in zGrid:
+        for y in yGrid:
+            for x in xGrid:
+                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 velocity float %d\n' %(1))
+    outfile.writelines('LOOKUP_TABLE default\n')
+
+    # write velocity
+    print("Writing velocity values to VTK file...")
+    for velocity in vel:
+        outfile.writelines('%10f\n' %velocity)
+
+    outfile.close()
+    print("Wrote velocity grid for %d points to file: %s" %(nPoints, outputfile))
+    return
+
+def rays2VTK(fnin, fdirout = './vtk_files/', nthPoint = 50):
+    '''
+    Writes VTK file(s) for FMTOMO rays from rays.dat
+
+    :param: nthPoint, plot every nth point of the ray
+    :type: integer
+    '''
+    def getDistance(angle):
+        PI = np.pi
+        R = 6371.
+        distance = angle / 180 * (PI * R)
+        return distance
+    
+    infile = open(fnin, 'r')
+    R = 6371
+    rays = {}
+    raynumber = 0
+    nPoints = 0
+
+    ### NOTE: rays.dat seems to be in km and radians
+
+    while True:
+        raynumber += 1 
+        firstline = infile.readline()
+        if firstline == '': break # break at EOF
+        shotnumber = int(firstline.split()[1])
+        nRayPoints = int(infile.readline().split()[0])
+        if not shotnumber in rays.keys():
+            rays[shotnumber] = {}
+        rays[shotnumber][raynumber] = []
+        for index in range(nRayPoints):
+            if index%nthPoint is 0 or index == nRayPoints:
+                rad, lat, lon = infile.readline().split()
+                rays[shotnumber][raynumber].append([getDistance(np.rad2deg(float(lon))), getDistance(np.rad2deg(float(lat))), float(rad) - R])
+            else:
+                dummy = infile.readline()
+
+    infile.close()
+
+    for shotnumber in rays.keys():
+        fnameout = fdirout + 'rays%03d.vtk'%(shotnumber)
+        outfile = open(fnameout, 'w')
+
+        nPoints = 0
+        for raynumber in rays[shotnumber]:
+            for ray in rays[shotnumber][raynumber]:
+                nPoints += 1
+
+        # write header
+        #print("Writing header for VTK file...")
+        print("Writing shot %d to file %s" %(shotnumber, fnameout))
+        outfile.writelines('# vtk DataFile Version 3.1\n')
+        outfile.writelines('FMTOMO rays\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 raynumber in rays[shotnumber].keys():
+            for raypoint in rays[shotnumber][raynumber]:
+                outfile.writelines('%10f %10f %10f \n' %(raypoint[0], raypoint[1], raypoint[2]))
+
+        outfile.writelines('LINES %15d %15d\n' %(len(rays[shotnumber]), len(rays[shotnumber]) + nPoints))
+
+        # write indices
+        #print("Writing indices to VTK file...")
+
+        count = 0
+        for raynumber in rays[shotnumber].keys():
+            outfile.writelines('%d ' %(len(rays[shotnumber][raynumber])))
+            for index in range(len(rays[shotnumber][raynumber])):
+                outfile.writelines('%d ' %(count))
+                count += 1
+            outfile.writelines('\n')
+
+    # outfile.writelines('POINT_DATA %15d\n' %(nPoints))
+    # outfile.writelines('SCALARS rays float %d\n' %(1))
+    # outfile.writelines('LOOKUP_TABLE default\n')
+
+    # # write velocity
+    # print("Writing velocity values to VTK file...")
+    # for velocity in vel:
+    #     outfile.writelines('%10f\n' %velocity)
+
+    # outfile.close()
+    # print("Wrote velocity grid for %d points to file: %s" %(nPoints, outputfile))
+

pylot/core/active/picking_script.py

@@ -0,0 +1,103 @@
+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 = "../rockeskyll_200615_270615/"
+else:
+    receiverfile = "Geophone_interpoliert_GZB" 
+    sourcefile = "Schusspunkte_GZB"
+    obsdir = "../GZB_26_06_15_01/"
+
+# 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)
+
+print '\n--- Finished script ---'
+print 'Elapsed time:', datetime.now()-starttime

pylot/core/active/seismicArrayPreperation.py

@@ -0,0 +1,665 @@
+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)

pylot/core/active/surveyPlotTools.py

@@ -0,0 +1,250 @@
+import matplotlib.pyplot as plt
+import math
+#from selectRegions import regions
+
+plt.interactive(True)
+
+def plotAllPicks(shot_dict, dist_med = None):
+    '''
+    Plots all picks over the distance between source and receiver. Returns (ax, region)
+    '''
+    dist = []
+    pick = []
+    snrloglist = []
+    for shot in shot_dict.values():
+        for traceID in shot.getTraceIDlist():
+            if shot.getPick(traceID) is not None: 
+                dist.append(shot.getDistance(traceID))
+                pick.append(shot.getPick(traceID))
+                snrloglist.append(math.log10(shot.getSNR(traceID)[0]))
+
+    ax = createPlot(dist, pick, snrloglist, label = 'log10(SNR)')
+    region = regions(ax, shot_dict)
+    if dist_med is not None:
+        ax = addDistMed(ax, dist_med)
+    ax.legend()
+
+    return ax, region
+
+def plotAllPicks_withCutOutTraces(shot_dict, dist_med = None):
+    '''
+    Plots all picks over the distance between source and receiver. Returns (ax, region)
+    '''
+    dist = []
+    pick = []
+    snrloglist = []
+    for shot in shot_dict.values():
+        for traceID in shot.getTraceIDlist():
+            if shot.getSNR(traceID)[0] > 3:
+                dist.append(shot.getDistance(traceID))
+                pick.append(shot.getPick_backup(traceID))
+                snrloglist.append(math.log10(shot.getSNR(traceID)[0]))
+
+    ax = createPlot(dist, pick, snrloglist, label = 'log10(SNR)')
+    region = regions(ax, shot_dict)
+    if dist_med is not None:
+        ax = addDistMed(ax, dist_med)
+    ax.legend()
+
+    return ax, region
+
+def createPlot(dist, pick, inkByVal, label):
+    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
+
+def addDistMed(ax, dist_med):
+    '''
+    Add shot dictionary containing the Median for several distancebins to the ax.
+    '''
+    x = []
+    y = []
+    for dist in dist_med.keys():
+        x.append(dist)
+        y.append(dist_med[dist])
+        
+    xy = sorted(zip(x,y))
+    x1 = [x for (x,y) in xy]
+    y1 = [y for (x,y) in xy]
+    
+    ax.plot(x1, y1, 'k', label = 'Median')
+
+    return ax
+
+class regions(object):
+    def __init__(self, ax, shot_dict):
+        self.ax = ax
+        self.shot_dict = shot_dict
+        self._x0 = []
+        self._y0 = []
+        self._x1 = []
+        self._y1 = []
+        self.shots_found = {}
+        self.shots_for_deletion = {}
+
+    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 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
+
+    def getShotDict(self):
+        return self.shot_dict
+
+    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... "
+
+        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]:
+                            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
+
+    def highlightPick(self, shot, traceID, annotations = True):
+        self.ax.scatter(shot.getDistance(traceID), shot.getPick(traceID), s = 50, marker = 'o', facecolors = 'none', edgecolors = 'm', alpha = 1)
+        if annotations == True:
+            self.ax.annotate(s = 's%s|t%s'%(shot.getShotnumber(), traceID), xy = (shot.getDistance(traceID), shot.getPick(traceID)), fontsize = 'xx-small')
+        self.ax.set_ylim(shot.getCut())
+
+    def plotTracesInRegion(self):
+        import matplotlib.pyplot as plt
+        count = 0
+        maxfigures = 20
+        # if len(self.shots_found) == 0:
+        self.findTracesInShotDict()
+
+        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 plotTracesInRegion_withCutOutTraces(self):
+        import matplotlib.pyplot as plt
+        count = 0
+        maxfigures = 20
+        # if len(self.shots_found) == 0:
+        self.findTracesInShotDict(picks = 'includeCutOut')
+
+        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 setCurrentRegionsForDeletion(self):
+        # if len(self.shots_found) == 0:
+        self.findTracesInShotDict()
+
+        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'
+
+    def markAllRegions(self, color = 'grey'):
+        from matplotlib.patches import Rectangle
+
+        for index in range(len(self._getx0())):
+            x0 = self._getx0()[index]
+            y0 = self._gety0()[index]
+            x1 = self._getx1()[index]
+            y1 = self._gety1()[index]
+
+            self.ax.add_patch(Rectangle((x0, y0), (x1 - x0), (y1 - y0), alpha=0.5, facecolor = color))
+            self.refreshFigure()
+
+    def markCurrentRegion(self, x0, x1, y0, y1, color = 'grey'):
+        from matplotlib.patches import Rectangle
+
+        self.ax.add_patch(Rectangle((x0, y0), (x1 - x0), (y1 - y0), alpha=0.1, facecolor = color))
+        self.refreshFigure()
+
+    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.refreshFigure()
+
+    def refreshFigure(self):
+        import matplotlib.pyplot as plt
+        plt.draw()