From b49206407ac3e10b964189f8f4067225838e5255 Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Sun, 18 Oct 2015 21:23:09 +0200 Subject: [PATCH] Merge using remote to resolve conflicts --- pylot/core/active/seismicshot.py | 107 ++++++++++++----------- pylot/core/active/surveyUtils.py | 143 +++++++++++++++++++++++++++++-- 2 files changed, 194 insertions(+), 56 deletions(-) diff --git a/pylot/core/active/seismicshot.py b/pylot/core/active/seismicshot.py index ecb602c2..df623479 100644 --- a/pylot/core/active/seismicshot.py +++ b/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), - nfac, (tnoise, tgap, tsignal), + (self.earliest[traceID], self.latest[traceID], tmp) = earllatepicker(self.getSingleStream(traceID), + 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)], - [min(stream[0].data), + ax1.plot([self.getPick(traceID), self.getPick(traceID)], + [min(stream[0].data), max(stream[0].data)], 'r', label = 'mostlikely') plt.legend() ax2 = plt.subplot(2,1,2, sharex = ax1) ax2.plot(hoscf.getTimeArray(), hoscf.getCF(), 'b', label = 'HOS') ax2.plot(hoscf.getTimeArray(), aiccf.getCF(), 'g', label = 'AIC') - ax2.plot([self.getPick(traceID), self.getPick(traceID)], - [min(np.minimum(hoscf.getCF(), aiccf.getCF())), + ax2.plot([self.getPick(traceID), self.getPick(traceID)], + [min(np.minimum(hoscf.getCF(), aiccf.getCF())), max(np.maximum(hoscf.getCF(), aiccf.getCF()))], 'r', label = 'mostlikely') ax2.plot([0, self.getPick(traceID)], @@ -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: diff --git a/pylot/core/active/surveyUtils.py b/pylot/core/active/surveyUtils.py index 909fed13..9fedd432 100644 --- a/pylot/core/active/surveyUtils.py +++ b/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 = value.replace(" ", "") + + value = parameterConfig.get('vars', parameter).split('\t')[0] return value +def setArtificialPick(shot_dict, traceID, pick): + for shot in shot_dict.values(): + shot.setPick(traceID, pick) + shot.setPickwindow(traceID, shot.getCut()) + def fitSNR4dist(shot_dict, shiftdist = 5): + import numpy as np dists = [] picks = [] snrs = [] @@ -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 + + +