#!/usr/bin/env python # # Script to create a bunch of nmo and reduction velocity horizons # for use in The Kingdom Suite. Aids in picking reflected and refracted # arrivals. Will try to produce a geoquest output file containing all # horizons (though these appear to need x and y coords, will need further testing). # Result - Works poorly # # # 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 import os import string as st import Numeric as N from optparse import OptionParser # from LatLongUTMconversion import LLtoUTM def getCmdLineArgs(): parser = OptionParser() parser.usage = """segyinfo2horizon.py -i input | -s input | -c input [options] [-o output] Script used to convert the segyinfo dump (or shottable) to an NMO or LMO horizon that can be imported into The Kingdom Suite. This horizon is useful for flattening the OBS data, and making it much easier to identify reflected (NMO) and refracted (LMO) arrivals. Currently individual files are created for each reduction velocity and NMO correction. These 'horizons' are also written to a single 'geoquest card image 7' file for quick import into the Kingdom Suite. The horizon files created when using the segy info dump file appear to work better, but must be imported individually as 'line trace time' format. These work better because they do not include missing traces, as the shottable does (use of shottable not recommended). Here is an example of how to create files for import into TKS, using the segyInfo dump from SeisWide (tempHeader). This command will output several different LMO and NMO curves based on what is seen in curveinfo.txt. Edit curveinfo.txt to create LMO and NMO curves with different velocities and depths. These files will be called 'obs_d_vel_red_1700.dat' for LMO of 1700 m/s for example. ./segyinfo2horizon.py -i tempHeader.txt -o obs_d.dat """ parser.add_option("-i", "--input", dest="inFile", \ help="Specify input file segyinfo dump format. Either this or the '-s' or '-c' " + \ "option must be provided [].", \ type="string", metavar="INPUT", default=None) parser.add_option("-s", "--shottable-input", dest="inShottable", \ help="Specify input file in shottable format [].", \ type="string", metavar="INPUT", default=None) parser.add_option("-c", "--claritas-input", dest="claritasInput", \ help="Specify input file in Claritas 'trprint' format. The only necessary " + \ "header field is 'offset'. [].", \ type="string", metavar="INPUT", default=None) parser.add_option("-l", "--line-name", dest="lineName", \ help="Specify line name used in output file [DefaultLine].", \ type="string", metavar="NAME", default="DefaultLine") parser.add_option("-r", "--curve-info", dest="curveinfo", \ help="Specify file with 'velocity, depth' info. Each line of the file " + \ "should have only the velocity in m/s and the depth in m, all whitespace " + \ "delimited. If the file does not exist, the file 'cruveinfo.txt' will be " + \ "created and used with default values [curveinfo.txt].", \ type="string", metavar="NAME", default="curveinfo.txt") parser.add_option("-o", "--output", dest="outFile", \ help="Specify output file in geoquest format [horizon_export.dat].", \ type="string", metavar="OUTPUT", default="horizon_export.dat") return parser.parse_args() def importShottable(inFile): tracesAndOffsets = [] traceCount = 0 for line in inFile.readlines(): fields = line.split() # only taking lines that can all be mapped as numbers try: fields = map(float, fields) except: continue if len(fields) == 0: continue # making a list of trace number, offset, lat, long tracesAndOffsets.append([traceCount + 1, fields[7], fields[5], fields[6]]) traceCount += 1 return tracesAndOffsets def importSegyinfo(inFile): tracesAndOffsets = [] for line in inFile.readlines(): # splitting on whitespace fields = line.split() if len(fields) == 13 : try: # Skipping lines that have an offset very close to 0 m if round(float(fields[4])*1000) == 0: continue tracesAndOffsets.append([int(fields[0]), float(fields[4]), 0, 0]) except ValueError: # Ignore the fact that fields[0] or float(fields[4]) doesn't work. pass return tracesAndOffsets def importTrprint(inFile): """import the Claritas trprint format file""" tracesAndOffsets = [] traceNo = 0 for line in inFile.readlines(): # skipping blank lines if not line: continue line = line.replace("="," ") # breaking up lines into a list of space delimited text fields = line.split() if not fields[0] == "%Trheader": continue # processing lines of the trprint output. 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]) tracesAndOffsets.append([traceNo + 1, offset / 1000.0, lat, longitude]) traceNo += 1 except: continue return tracesAndOffsets def makeLmoCurves(lmoVelocity, offsetArray): lmoArray = N.copy.copy(offsetArray) # reduction velocity is in m/s, and offset in km. Converting to km/s lmoArray[0:,1] = N.absolute(lmoArray[0:,1]) / (lmoVelocity / 1000) return lmoArray def makeNmoCurves(nmoVel, nmoDepth, offsetArray): # Calculate the NMO correction. Assumes 1 way time to OBS in water. NmoArray = N.copy.copy(offsetArray) offset = offsetArray[0:,1] * 1000 # new array is 'shot, time' NmoArray[0:,1] = (((nmoDepth**2 / nmoVel**2) + (offset**2 / nmoVel**2))**0.5 - nmoDepth / nmoVel) return NmoArray def writeXyFile(dataArray, filename, lineName): # writing out the data to a file with a decent format. This is easily imported by The Kingdom Suite. outputObj = open(fileName, "w") for point in range(0, N.shape(dataArray)[0]): outputObj.write("%s %05d %6.3f\n" % (lineName, dataArray[point,0], dataArray[point,1])) outputObj.close() def writeGeoquest(inArray, fileName, lineName, velList, depthList, horizType): # writing a Geoquest Card Image 7 format file. TKS appears to have problems if there are no # coordinates in this file. horizonIdx = 0 for horizon in inArray: # printing the header(s) for each horizon if horizType == "nmo": horizonString = "nmo_v" + str(velList[horizonIdx]) + "_z" + str(depthList[horizonIdx]) else: horizonString = "velocity_reduct_" + str(velList[horizonIdx]) outputObj.write("PROFILE %s TYPE 1 3 \n" % (horizonString)) outputObj.write("SNAPPING PARAMETERS 2 200 2\n") for point in horizon: strTime = ("%7.2f" % (point[1] * 1000)) # time is in ms. strShotNo = ("%4.2f" % point[0]) strShotNoInt = str(int(point[0])) strAmplitude = ("%08.5f" % 0.0) shotX = point[2] shotY = point[3] # rjust requires strings, no ints or floats. First two zeros correspond to coordinates - none for now. outputObj.write(" %1.8E %1.8E %s %s %s %s 2 %s\n" % \ (shotX, shotY, st.rjust(strTime, 14), st.rjust(strShotNo, 9), st.rjust(strShotNoInt, 14), st.rjust(strAmplitude, 13), lineName)) # printing footers for each horizon outputObj.write("EOD %s\n" % horizonString) horizonIdx += 1 if __name__ == "__main__": # getting the command line options from the parser (options,args)= getCmdLineArgs() if not (options.inFile or options.inShottable or options.claritasInput): print "Error: please specifiy an input file. 'Type segyinfo2horizon.py -h' for detailed help." #print options, args sys.exit() # importing the input data, each function returns a python list if options.inFile: inputObj = open(options.inFile) traceAndOffset = importSegyinfo(inputObj) elif options.inShottable: inputObj = open(options.inShottable) traceAndOffset = importShottable(inputObj) elif options.claritasInput: inputObj = open(options.claritasInput) traceAndOffset = importTrprint(inputObj) # converting the list no a numeric array traceAndOffset = N.array(traceAndOffset) # reading values of lmo and nmo correction from file (s) if it exists, if not then make it. if not os.path.isfile(options.curveinfo): print "curve info file does not exist, creating file", options.curveinfo, "with default values" outputObj = open(options.curveinfo, "w") outputObj.write("# this file automatically created by segyinfo2horizon.py.\n") outputObj.write("# Please customize values below. Each line corresponds to a single\n") outputObj.write("# curve that will be created. Single numbers on a line are treated\n") outputObj.write("# as LMO reduction velocities, while two numbers are treated as\n") outputObj.write("# NMO velocity and depth. (all in m and m/s)\n") outputObj.write("# LMO velocities (m):\n") outputObj.write("%d\n" % (1550)) outputObj.write("# NMO velocities (m/s) and depths (m):\n") outputObj.write("%d %d\n" % (1500, 1800)) outputObj.write("%d %d\n" % (1700, 2000)) outputObj.close() lmoList = [] velList = [] depthList = [] curveInfoFile = open(options.curveinfo, "r") for velDepthLine in curveInfoFile.readlines(): # if something has a '#' mark in the line, we ignore the rest of the line nonComments = velDepthLine.split('#')[0] fields = nonComments.split() # skipping blank lines if len(fields) == 0: continue # single float/int value: treat as LMO reduction. if len(fields) == 1: try: fields = map(float, fields) lmoList.append(fields[0]) except: pass # two floats/ints, treat as NMO reduction. elif len(fields) == 2: try: fields = map(float, fields) velList.append(fields[0]) depthList.append(fields[1]) except: pass # Calculating and writing out the velocity reduction files. lmoArray = N.array([]) traceCount = None for velocity in lmoList: # using a seperate function to calc the traveltime reductions. lmoHorizon = makeLmoCurves(velocity, traceAndOffset) traceCount = len(lmoHorizon) # write out the lmoArray to a file, with descriptive filename fileName = os.path.splitext(options.outFile)[0] + "_vel_red_" + \ str(int(velocity)) + os.path.splitext(options.outFile)[1] # making a large array for writing to the geoquest format later on. if not lmoArray: lmoArray = lmoHorizon else: lmoArray = N.concatenate([lmoArray, lmoHorizon]) writeXyFile(lmoHorizon, fileName, options.lineName) # reshaping array so the smallest dimension is 2 (trace, time), followed by n times the number of traces, # followed by the number of horizons (LMO velocities/curves). lmoArray = N.reshape(lmoArray, (len(lmoList), traceCount, 4)) # # opening geoquest file object - not doing anymore, these curves are imported like crap. # nmoFileName = os.path.splitext(options.outFile)[0] + "_geoquest" + os.path.splitext(options.outFile)[1] # outputObj = open(nmoFileName, "w") nmoArray = N.array([]) nmoIdx = 0 traceCount = None # creating velocity and depth lists while nmoIdx < len(velList): velocity = velList[nmoIdx] depth = depthList[nmoIdx] # velList.append(velocity) # depthList.append(depth) nmoHorizon = makeNmoCurves(velocity, depth, traceAndOffset) traceCount = len(nmoHorizon) # making a large array for writing to the geoquest format later on. if not nmoArray: nmoArray = nmoHorizon else: nmoArray = N.concatenate([nmoArray, nmoHorizon]) # write out the nmoArray to an xy file, with descriptive filename fileName = os.path.splitext(options.outFile)[0] + "_nmo_vel_" + str(int(velocity)) + \ "_depth_" + str(int(depth)) + os.path.splitext(options.outFile)[1] writeXyFile(nmoHorizon, fileName, options.lineName) nmoIdx += 1 # reshaping the NMO array to be the same as LMO array. nmoArray = N.reshape(nmoArray, (len(velList), traceCount, 4)) # only the shottable has input coords, so ignoring all others. if options.inShottable: # converting the lat / long coords to UTM X and Y. Function not vectorized, so doing it the SLOWWW way for profile in lmoArray: for point in profile: (utmZone, point[2], point[3]) = LLtoUTM(23, point[2], point[3]) for profile in nmoArray: for point in profile: (utmZone, point[2], point[3]) = LLtoUTM(23, point[2], point[3]) # # writing out the geoquest format file for both velocity reduction and NMO horizons. # # Not doing this anymore, TKS imports it like crap # writeGeoquest(lmoArray, nmoFileName, options.lineName, lmoList, depthList, "lmo") # writeGeoquest(nmoArray, nmoFileName, options.lineName, velList, depthList, "nmo") # # outputObj.close()