#!/usr/bin/env python # This is a program that converts picks from Geoquest7 to Zelt's tx.in format. # It is intended for OBS data # # Originally the trace in the exported file would correspond to the shot # in the shottable. However, Kingdom suite is renumbers the traces sequentially # therefore if shots are missing, things get screwed up. Instead, the proper # relationship of trace to shot is derived from the "dump more trace info" option # of "Segy Info" in Seiswide. This makes this program very narrow in its scope. # # The Geoquest image format was used because it exports lots of useful information, and # can handle multiple layers. Therefore only one conversion is needed for all horizons. # # I have decided to not use any traces with exactly 0 m offset, because these are probably # garbage. Dummy traces added to the SEGY are treated as having 0 m offset, therefore # messing up the analysis. TODO: clean it up so this logic is not repeated for each output # filetype. # # 22 march 2005. Added logic to the 'transformations' to adjust the by X_OFFSET. # # # This program is free software; you can redistribute it and/or # modify it under the terms of the GNU program General Public # License as published by the Free Software Foundation; either # version 2 of the License, or (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # program General Public License for more details. # # You should have received a copy of the GNU program General Public # License along with this program; if not, write to the Free Software # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA # # $Id: geoquest2zelt.py,v 1.34 2004/04/12 20:25:07 chris Exp $ import sys, getopt, fileinput, unittest import StringIO from optparse import OptionParser usage = """ An all platform program intended for conversion from the 'geoquest CARD IMAGE 7' export (kingdom Suite default) to the Zelt format. You will require both the exported file from Kingdom Suite, and the 'dump more trace info' file from SeisWide (.tempheader.txt). The only required option is the input pick file, all others are optional. This program is distributed under the GNU General Public License. Synopsis:'geoquest2zelt.py [OPTIONS]... -i INPUT.TXT' Examples:Converting 'foo.dat' to tx.in, with the Segy info file '.tempHeader.txt'. In DOS, leave out the './': './geoquest2zelt.py -i foo.dat' Converting foo.dat to zelt.in, with pick uncertainty of 40 ms, removing a 6 km/s reduction velocity (therefore inverse reduction), and applying a static shift of 0.5 s: './geoquest2zelt -e 0.4 -I -r 6 -t 0.5 --input=foo.dat \ --output=zelt.in' Using the Claritas trprint file 'trprint.txt' instead of a segyInfo file. This will output to tx.in: './geoquest2zelt.py -i foo.dat -c trprint.txt' NOTE: Take extreme caution if you have used processes such as 'missing' which vary the number of traces in a section. This program renumbers the offsets to increase monotonically from one, just as TKS does. If the number of traces in the segy file doesn't match what is in the segyinfo or claritasinfo file, this will certainly give incorrect values. The best plan is to write out the trprint file in the same job the filtered segy file is written in (prior to importing into TKS). That way the files should always be in sync. """ def getCmdLineArgs(args): parser = OptionParser() parser.add_option("-i", "--input", dest="geoFile", \ help="Specify input horizon file Must be in Geoquest Card Image 7 " + \ "format, the Kingdom Suite default. Horizons must start with the " + \ "shallowest first. Mandatory input.", \ type="string", metavar="INPUT", default=None) parser.add_option("-o", "--output", dest="outFile", \ help="Output file in 'tx.in' format, default is 'tx.in'.", \ type="string", metavar="OUTPUT", default="tx.in") parser.add_option("-s", "--segy-info", dest="segyInfo", \ help="'segy info' file output from SeisWide. Used to match" + \ "'shot numbers' with 'trace numbers' and 'offset'. Can be any " + \ "whitespace delimited text file containing 'trace', 'shot' (FFID), " + \ "and 'offset' (km) as columns 1, 2, and 5 respectively (no other " + \ "columns used). Trace number is always monotonically renumbered from 1," + \ "to match the behaviour of TKS. If this " + \ "is not specified, the '-c' option must be.", \ type="string", metavar="SEGY_INFO", default=None) parser.add_option("-c", "--claritas-info", dest="claritasInfo", \ help="Output from Claritas 'trprint'. Used to match shot numbers " + \ "with trace numbers and offset. Must have 'offset' header output in file in metres.", \ type="string", metavar="CLARITAS_INFO", default=None) parser.add_option("-x", "--x-offset", dest="obsSourceOffset", \ help="X offset of OBS from origin of Zelt velocity model. " + \ "Default is 0 km.", \ type="float", metavar="X_OFFSET", default=0.0) parser.add_option("-e", "--error", dest="pickError", \ help="Uncertainty of picks in 'pickfile', default is 0.010 s.", \ type="float", metavar="PICK_ERROR", default=0.01) parser.add_option("-t", "--time-shift", dest="timeshift", \ help="Static shift to apply to travel times in seconds. " +\ "Default is 0 s.", \ type="float", metavar="TIME_SHIFT", default=0.0) parser.add_option("-r", "--reduction", dest="reductVelocity", \ help="Apply a forward (optionally inverse) reduction velocity to the picks, " + \ "in km/s. Must be a positive number. Default is 0 km/s for no reduction.", \ type="float", metavar="REDUCTION_VELOCITY", default=0.0) parser.add_option("-I", "--inverse", action="store_true", dest="inverse", \ help="Apply the inverse of the reductin velocity specified in '--reduction'. " + \ "This is a flag, so no arguments.", \ metavar="INVERSE_FACTOR") parser.add_option("-X", "--xyformat", action="store_true", dest="xyformat", \ help="Write out a SeisWide XY format 'offset, traveltime, phase' instead of "+ \ "the default Zelt format. Useful for importing picks into Seiswide. No arguments.", \ metavar="XY_FORMAT") parser.add_option("-l", "--obsloc", action="store_true", dest="obslocFormat", \ help="Write out an ascii file of 'shotNumber, traveltime' useful for " + \ "many things, especially the OBS relocation program ObsLoc. No arguments.", \ metavar="OBSLOC_FORMAT") return parser.parse_args(args) class FileReader: def __init__(self, file): if type(file) == str: self.fd = open(file, "rb") else: # Assume it is a file self.fd = file def readFile(self): for line in self.fd.readlines(): # There is probably a better way to strip newlines line = line.replace("\n", "") line = line.replace("\r", "") if len(line) != 0: self.processLine(line) def processLine(self, line): pass class GeoquestReader(FileReader): def __init__(self, file, phases): FileReader.__init__(self, file) self.phases = phases self.phaseNum = 0 def processLine(self, line): line = line.split() if line[0] == "EOD" : self.phases.append(Phase()) self.phaseNum += 1 return None else: if (not line) or (len(line) < 8): return None lineDict = {"travelTime" : float(line[2]), "traceNum" : int(line[4]), "phaseNum" : self.phaseNum} if len(self.phases) <= self.phaseNum: self.phases.append(Phase()) self.phases[self.phaseNum].append(lineDict) class SegyInfoReader(FileReader): def __init__(self, file, phases): FileReader.__init__(self, file) self.phases = phases # using a counter to renumber trace number instead of # using SeisWide's 'trace number' which seems to change # between starting from 0 and starting from 1. This new # behaviour matches what TKS does. self.traceCount = 0 def processLine(self, line): try: fields = map(float, line.split()) if len(fields) != 13: return None #traceNum = int(fields[0]) # old 'seiswide' trace number traceNum = self.traceCount + 1 # new renumbered version shotNum = int(fields[1]) offset = float(fields[4]) updateLines = Phase.getMatchingLines(self.phases[0], traceNum) if updateLines: for updateLine in updateLines: updateLine["offset"] = offset updateLine["shotNum"] = shotNum self.traceCount += 1 except ValueError: # Not converted correctly, ignore return None class ClaritasInfoReader(FileReader): def __init__(self, file, phases): FileReader.__init__(self, file) self.phases = phases # using a counter to renumber trace number instead of # using SeisWide's 'trace number' which seems to change # between starting from 0 and starting from 1. This new # behaviour matches what TKS does. self.traceCount = 0 def processLine(self, line): try: line = line.replace("="," ") # breaking up lines into a list of space delimited text fields = line.split() if fields[0] != "%Trheader": return None # processing lines of the trprint output. traceNum = self.traceCount + 1 offset = 0.0 lat = 0.0 longitude = 0.0 # only taking lines that can be mapped as floating point values (no text) try: for fieldNo in range(0, len(fields)): if fields[fieldNo] == "OFFSET": offset = float(fields[fieldNo + 1]) / 1000 updateLines = Phase.getMatchingLines(self.phases[0], traceNum) if updateLines: for updateLine in updateLines: updateLine["offset"] = offset #updateLine["shotNum"] = shotNum self.traceCount += 1 except: return None except: # Not converted correctly, ignore return None class OutputFile: """Output to an XY format file. override writePhases to output other formats""" def __init__(self, options): self.options = options # making an object from the string name of the output file self.outFile = open(self.options.outFile, "w") def close(self): self.outFile.close() def writePhases(self, phases): phaseCount = 0 for phase in phases: for line in phase: #Discarding any offsets of zero offset if line["offset"] == options.obsSourceOffset: continue #Discarding any null traveltimes, hardcoded. if line["travelTime"] == -0.999: continue try: self.outFile.write(" %6.3f %6.3f %d curve%d\n" %\ (line["offset"], line["travelTime"], phaseCount + 1, phaseCount + 1)) except KeyError: print ("Warning: no offset for trace no %d" % line["traceNum"]) phaseCount += 1 self.outFile.flush() class ObslocOutputFile(OutputFile): def writePhases(self, phases): phaseCount = 0 xyList = [] for i in range(len(phases)): xyList.append([]) for phase in phases: for line in phase: #Discarding any offsets of zero offset if line["offset"] == options.obsSourceOffset: continue #Discarding any null traveltimes, hardcoded. if line["travelTime"] == -0.999: continue # Todo: sort this data from smallest to largest shot number - needed by obsloc xyList[phaseCount].append(line) # # outputting shotNumber (recordnum) instead of tracenumber for ObsLoc. # try: # self.outFile.write("%d %.2f\n" %(line["shotNum"], line["travelTime"] * 1000.0)) # except KeyError: # print ("Warning: no offset for trace no %d" % # line["traceNum"]) xyList[phaseCount].sort(cmpTraceFunc) phaseCount += 1 # Writing out the sorted xy list. Being dumb, and not handling exceptions. for phase in xyList: for line in phase: self.outFile.write("%d %.2f\n" %(line["shotNum"], line["travelTime"] * 1000.0)) self.outFile.flush() # new function for sorting based on shot numbers def cmpTraceFunc(line1, line2): """Return 1 if line1 > line2, -1 if less, 0 if equal.""" difference = line1["shotNum"] - line2["shotNum"] if difference > 0: return 1 elif difference < 0: return -1 else: return 0 def cmpLineFunc(line1, line2): """Return 1 if offset of line1 > line2, -1 if less, 0 if equal.""" difference = line1["offset"] - line2["offset"] if difference > 0: return 1 elif difference < 0: return -1 else: return 0 class ZeltOutputFile(OutputFile): def writePhases(self, phases): phaseNum = 0 # Sort into negative and postive groups still ordered by phases negList = [] posList = [] for i in range(len(phases)): negList.append([]) posList.append([]) for phase in phases: for line in phase: try: #Discarding any offsets of zero offset if line["offset"] == options.obsSourceOffset: continue #Discarding any null traveltimes (-999 ms), hardcoded. if line["travelTime"] == -0.999: continue if line["offset"] < options.obsSourceOffset: negList[phaseNum].append(line) elif line["offset"] > options.obsSourceOffset: posList[phaseNum].append(line) except KeyError: print("Warning: no offset for trace no %d" % line["traceNum"]) negList[phaseNum].sort(cmpLineFunc) posList[phaseNum].sort(cmpLineFunc) phaseNum += 1 if len(negList[0]) > 0: #1st sentinel. Right justifying certain values with rjust self.outFile.write("%s %06.3f %05.3f %d\n" % (str("%6.3f" % self.options.obsSourceOffset).rjust(10), -1, 0, 0)) self.writeBlock(negList) if len(posList[0]) > 0: self.outFile.write("%s %6.3f %05.3f %d\n" % (str("%6.3f" % self.options.obsSourceOffset).rjust(10), 1, 0, 0)) self.writeBlock(posList) # End sentinel self.outFile.write(" %5.3f %05.3f %05.3f %d\n" % (0, 0, 0, -1)) def writeBlock(self, block): layer = 0 for phase in block: for line in phase: self.outFile.write(\ "%s %s %05.3f %s\n" % (str("%6.3f" % line["offset"]).rjust(10), str("%6.3f" % line["travelTime"]).rjust(5), self.options.pickError, str(layer + 1).rjust(9))) layer += 1 class Phase: """Basically is a list that keeps track of all lines matching a particular traceNum in a (class) attribute that maps traceNum to a list of lines. Does not update when a line is deleted. """ traceToLines = {} def __init__(self): self.lineList = [] def __len__(self): return len(self.lineList) def append(self, line): self.lineList.append(line) traceNum = line["traceNum"] self.traceToLines[traceNum] = self.traceToLines.get(traceNum, []) + [line] def __getitem__(self, index): return self.lineList[index] def __iter__(self): return iter(self.lineList) def __str__(self): return str(self.lineList) def sort(self, cmpFunc): return self.lineList.sort(cmpFunc) def getMatchingLines(self, traceNum): #two tabs to the following item: if self.traceToLines.has_key(traceNum): return self.traceToLines[traceNum] else: return None def applyTransforms(phases, options): for phase in phases: # Handling the case where there is a trace in the pick file, but no corresponding trace/offset # in the Segyinfo file. Assigning 0 km offset so the other logic doesn't freek out. A warning is # printed and the trace later discarded (because of 0 km offset). Stderr may be a better way to go, oh well. for line in phase: try: if line["offset"]: pass except: print ("Warning: no corresponding trace/shot/offset for trace no %d in trace info file." % line["traceNum"]) line["offset"] = 0.0 phase.sort(cmpLineFunc) for line in phase: travelTime = line["travelTime"] travelTime = travelTime / 1000 + \ options.timeshift # optional timeshift, default is 0. # I've learned that you have to add the offset to the picks, # in addition to the 'offset' tags. if options.obsSourceOffset: line["offset"] = line["offset"] + options.obsSourceOffset if options.reductVelocity != 0.0: # performing the reduction velocity correction reduction = abs(line["offset"] / options.reductVelocity) # handling the inverse case if options.inverse: reduction = reduction * -1 # applying the reduction travelTime = travelTime - reduction line["travelTime"] = travelTime if __name__ == "__main__": (options, args)= getCmdLineArgs(sys.argv[1:]) # Checking for mandatory command line arguments if options.geoFile == None: print("Error: you must supply an input file (--input). Usage follows\n") sys.exit(usage) # Our list of phases - with an initially empty 1st phase phases = [] geoFile = GeoquestReader(options.geoFile, phases) geoFile.readFile() # Add the offset if options.segyInfo: segyInfoFile = SegyInfoReader(options.segyInfo, phases) segyInfoFile.readFile() elif options.claritasInfo: claritasInfoFile = ClaritasInfoReader(options.claritasInfo, phases) claritasInfoFile.readFile() applyTransforms(phases, options) # for line in phases: # for point in line: # print point if options.xyformat: outFile = OutputFile(options) elif options.obslocFormat: outFile = ObslocOutputFile(options) else: outFile = ZeltOutputFile(options) outFile.writePhases(phases) outFile.close() # Printing a final statement, saying "we're done!": print ("\nZelt file successfully output to: %s\n" % (options.outFile))