#!/usr/bin/env python # A program to convert a velocity model from Rayinvr to # to a Claritas style .nmo text format or xyz text file. # # 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 import sys import fileinput from string import * from optparse import OptionParser def getCmdLineArgs(): parser = OptionParser() parser.usage = """rayinvr2claritas.py --input [input] [options] --output [output] A script to convert from a Rayinvr velocity model to a Claritas velocity model. Velocity(time) profiles for specified distances must have been calculated and exist in 'vm.out' (run the Rayivnr program 'vmodel', with the paramters 'ivp=1', 'izort=1', 'xvp=offset1,offset2,...,offsetN', see rayinvr docs for details). These velocity profiles will be converted to format that Claritas will read. The program will take the 'vm.in' format as input, or a GMT style XYZ file (offset km, traveltime s, velocity km/s). It can also write out either a Claritas format .nmo file (default), or a GMT style .xyz format ascii file. The input xyz file should be sorted starting with the lowest offset (or trace number), and starting at the top of the model (0 sec two way time) progressing to the bottom of the model, then onward to higher offsets or trace numbers. Note: Claritas only uses velocity files in terms of travel time. Beware that Claritas has problems using negative offsets in velocity models. A solution for this is to use the '--trace-file' option, providing a file that contains 'trace, offset' pairs. This program will then write out a file containing trace numbers instead of offsets. Note: this program discards offsets of exactly -1 km, because this will cause an error in Claritas. Using trace numbers is a way to avoid this. Examples: This simple example converts 'vm.out' to 'vm.nmo', keeping OFFSET as the primary key. rayinvr2claritas.py --input vm.out --output vm.nmo Here is a shortened version of this command (the default output is vm.nmo): rayinvr2claritas.py -i vm.out This example converts the xyz style input file 'vm.xyz' to an xyz output, but writes out 'trace, time, velocity' instead of 'offset, time, velocity'. rayinvr2claritas.py --input vm.xyz --xyz-input --trace-file horizonOffset3.txt --xyz-output --output vm2.xyz """ parser.add_option("-i", "--input", dest="inFile", \ help="Specify input file in the Rayinvr 'vm.out' format, or GMT " + \ "style .xyz format, see the -x option. Must contain velocity " + \ "profiles in velocity vs. time. This is a mandatory option", \ type="string", metavar="INPUT", default=None) parser.add_option("--xyz-input", action="store_true", dest="xyzInput", \ help="Specify that input file is in GMT style .xyz format, not " + \ "default 'vm.out' format.", \ metavar="XYZ_INPUT") parser.add_option("-s", "--shift", dest="dataShift", \ help="Shift to add to all input traveltimes (or depths) in seconds. [0.0]", \ type="float", metavar="SHIFT", default=0.0) parser.add_option("-d", "--depth", action="store_true", dest="inputIsDepth", \ help="Specify that input file is in 'offset(km), depth(km), velocity(km/s)', " + \ "not 'offset(km), two-way time(s), velocity(km/s)'.", \ metavar="DEPTH_INPUT") parser.add_option("--trace-file", dest="traceFile", \ help="Specify a 'trace offset' file, to write a file that uses trace numbers " + \ "instead of offsets. Input must be an xy file, whitespace delimited, with trace " + \ "number as the first column, and offset (km) as the second column. This will " + \ "pick the trace it sees with the closest offset to the velocity profile(s) offset.", \ type="string", metavar="TRACE_FILE", default=None) parser.add_option("-o", "--output", dest="outFile", \ help="Specify output file in the default Claritas .nmo format, or GMT " + \ "style .xyz format, see the -X option [vm.nmo].", \ type="string", metavar="OUTPUT", default="vm.nmo") parser.add_option("--xyz-output", action="store_true", dest="xyzOutput", \ help="Output file should be in GMT style .xyz format, not default Claritas " + \ " .nmo format.", metavar="XYZ_INPUT") return parser.parse_args() def readRayinvr(VelProfile): # load the vm.out file into a list. # Checking of input data, should be 5 or 7 columns of data. inputData = [] xoffset = -99999 for line in fileinput.input(VelProfile): fields = line.split() # Obtaining the x offset try: if str(fields[0]) == "x" and str(fields[1]) == "position:": xoffset = float(fields[2]) #print xoffset except: continue # Only taking lines with 5 or 7 columns if len(fields) < 5: continue if len(fields) > 7: continue if len(fields) == 6: continue # Eliminating all lines that dont have a number in the # second column try: float(fields[1]) except ValueError: continue # Deleting the "layer# x" fields if str(fields[0]) == "layer#": del fields[0:2] inputData.append([float(xoffset), float(fields[1]), float(fields[3])]) # rearranging data so it is in the .xyz format xyzData = [] for idx in range(1, len(inputData), 2): zOffset, z1, z2 = inputData[idx - 1] vOffset, vp1, vp2 = inputData[idx] xyzData.append([zOffset, z1, vp1]) xyzData.append([vOffset, z2, vp2]) return xyzData def nestXyz(xyzInput): #print xyzInput allProfiles = [] singleProfile = [] # assigning each offset to its seperate list inside a larger list (list in a list) for point in range(0, len(xyzInput)): currentOffset = xyzInput[point][0] prevOffset = xyzInput[point -1][0] if currentOffset == prevOffset: singleProfile.append(xyzInput[point]) # if its a different offset, start a new sub-list and put values in it else: singleProfile = [] singleProfile.append(xyzInput[point]) allProfiles.append(singleProfile) # if there was only one offset, then it should still be appended to the superlist # for simplicity later on. if len(allProfiles) == 0: allProfiles.append(singleProfile) return allProfiles def shiftXyz(xyzNested, dataShift): for profile in xyzNested: for sample in profile: # adding a time shift to all time values sample[1] = sample[1] + dataShift # never want the first sample in a profile to be non zero. for profile in xyzNested: profile[0][1] = 0.0 return xyzNested def depthToTime(xyzNested): timeList = [] for profile in xyzNested: # print profile # assumes model starts at 0 m depth, going downward (+). totalTime = 0.0 timeList.append(profile[0]) for point in range(1, len(profile)): currentDepth = float(profile[point][1]) previousDepth = float(profile[point-1][1]) currentVel = float(profile[point][2]) previousVel = float(profile[point-1][2]) # depth of the interval intervalDepth = currentDepth - previousDepth # since velocity is a linear gradient, we can use the average avgIntVelocity = (currentVel + previousVel) / 2 # conversion of depth to 2 way time twoWayTime = (intervalDepth / avgIntVelocity) * 2 totalTime = totalTime + twoWayTime # adding to the new list, which will replace the old one timeList.append([profile[point][0], totalTime, profile[point][2]]) return timeList def offsetToTrace(xyzNested, traceFile): traceList = [] for line in file(traceFile, "r").readlines(): fields = line.split() # discarding lines without 3 columns if len(fields) != 2: continue # discarding any headers or lines that can't be read as such try: trace = int(fields[0]) offset = float(fields[1]) except: continue traceList.append([trace, offset]) # using loops (slow) to go through all the data points, replacing offset with trace xyzNestedTraces = [] currentOffset = None currentTrace = None for profile in xyzNested: # for sample in profile: # print sample for idx in range(0, len(profile)): #print profile[idx] leastOffset = 10000 # calculating closest trace number for each new offset if profile[idx][0] != currentOffset: currentOffset = profile[idx][0] for lineNum in range(0, len(traceList)): if abs(profile[idx][0] - traceList[lineNum][1]) < leastOffset: leastOffset = abs(profile[idx][0] - traceList[lineNum][1]) currentTrace = traceList[lineNum][0] #print currentTrace, currentOffset xyzNestedTraces.append([currentTrace, profile[idx][1], profile[idx][2]]) return xyzNestedTraces def writeXyz(inputList, outFile, xyzHeader): outFileObj = open(outFile, "w") outFileObj.write("XYZ output from Rayinvr2claritas.py\n" + \ xyzHeader + "\n") for Profiles in inputList: for point in Profiles: outFileObj.write("%.2f %.2f %.2f\n" % (point[0], point[1], point[2])) outFileObj.close() # nice message at the end telling where the file was written. sys.stdout.write("xyz output written to: " + options.outFile + "\n") def writeNmo(inputList, outFile, nmoHeader): outFileObj = open(outFile, "w") # Write out the .nmo header fluff # if the first column of the input file doesn't represent offset, this should # be changed accordingly. outFileObj.write("NMO\n" + \ "PRIMARY KEY : " + nmoHeader + "\n" + \ "SECONDARY KEY : ???\n" + \ "INTERPOLATION KEY : Interpolate/End\n" + \ "Created:\n" + \ "|KEYVAL|{T1 |V1} |{T2 |V2} |{T3 |V3} |{T4 |V4}... |\n") # # Taking the average velocity of each layer. Perhaps unnecessary, so excluded. # for i in range(len(inputData) - 1): # inputData[i][1] = (inputData[i][1] + inputData[i + 1][1]) / 2 # For each clump of 4 lines in the input file, write them to one line # in the output file. for profile in inputList: inputLen = len(profile) step = 4 for i in range(0, inputLen, step): # Have to skip offsets of exactly "-1", because Claritas can't cope. #print int(profile[i][0]) if int(profile[i][0]) == -1: continue # For the first line only, write out the "cdp" (or offset), else write an equal # amount of space. if i == 0: offset = str(int(profile[i][0])) whiteSpace = abs(7 - len(offset)) # must stick in some blank space so everything is alligned vertically # Changed 'offset' from an int to a float, because Claritas can handle it. # changing everything to strings, so its left justified. outFileObj.write(" %s%s" %(offset, " " * whiteSpace)) else: # sticking in space to fill in the column in which offset isn't written. - crazy format. outFileObj.write(" ") for j in range(step): if i + j < inputLen: # Writing of the nmo format, columns have to be left justified, and therefore # have to be converted to strings, hence the ugly syntax. outFileObj.write("%s %s " %((str(int(profile[i + j][1]*1000))).ljust(6),\ (str(int(profile[i + j][2]*1000))).ljust(6))) outFileObj.write("\n") sys.stdout.write("Claritas nmo output written to: " + outFile + "\n") outFileObj.close() # Main part of the program if __name__ == "__main__": # getting the command line options from the parser (options,args)= getCmdLineArgs() # if no input specified, give helpful error if not options.inFile: print "Error: must specify input file.\n" + \ "type 'rayinvr2claritas --help' for help" sys.exit() # reading the vm.out file (if its not an xyz format) inputData = [] if not options.xyzInput: inputData = readRayinvr(options.inFile) # if its an .xyz input file, load that into the same list if options.xyzInput: inputData = [] # Checking of input data, should be 3 columns of numbers (offset, depth[or tt], velocity), for line in fileinput.input(options.inFile): fields = line.split() if len(fields) != 3: continue try: fields = map(float, (fields)) except ValueError: continue inputData.append(fields) # Now that either input is in the same format, convert from a single list to a nested list. xyzNested = nestXyz(inputData) # perform any transformation on the data. Anything new is easy to add in the above method. if options.dataShift: xyzNested = shiftXyz(xyzNested, options.dataShift) # if data is in depth, convert to two way time if options.inputIsDepth: xyzNested = depthToTime(xyzNested) # nesting again since it returned a single list xyzNested = nestXyz(xyzNested) # default 'headers' for xyz and nmo files. xyzHeader = "[x offset(km), 2wayTime (s), velocity(km/s)]" nmoHeader = "OFFSET" # if offsets are to be converted to trace numbers if options.traceFile: # converting to trace numbers xyzNested = offsetToTrace(xyzNested, options.traceFile) # nesting the new trace list by trace number xyzNested = nestXyz(xyzNested) # changing the headers on either output file, to reflect that its now traces. xyzHeader = "[trace number, 2wayTime (s), velocity(km/s)]" nmoHeader = "LINE" # if an xyz output file is desired if options.xyzOutput: writeXyz(xyzNested, options.outFile, xyzHeader) # if the default .nmo output file is desired if not options.xyzOutput: writeNmo(xyzNested, options.outFile, nmoHeader)