convert_SU_to_ASCII.py

#!/usr/bin/env python
#
# reads in an SU-format file generated by SPECFEM and outputs ASCII files for each trace
#
"""
usage: ./convert_SU_to_ASCII.py file_in

  with file_in - input SU-format file, e.g., OUTPUT_FILES/0_dz_SU

"""
from __future__ import print_function

import sys,os
import numpy as np

# based on obspy reader
try:
    from obspy import read
except:
    print("Failed importing obspy. Please make sure to have version >= 1.1.0 working properly.")
    sys.exit(1)

# arguments
try:
    file_in = sys.argv[1]
except:
    print(__doc__)
    sys.exit(1)

# for example: ./convert_SU_to_ASCII.py OUTPUT_FILES/0_dz_SU
#              then filename = './OUTPUT_FILES/0_dz_SU' and single traces
#              will be outputted as 0_dz_SU.ascii_0, 0_dz_SU.ascii_1, etc.
out_file = file_in + '.ascii'

# by default: uses Seismic Unix format to read in file.
#             one could also use SEG-Y files and change this setting.
read_segy = False
read_su = True

# endianness:
# On the Macintosh platform, PowerPC-based Macintosh computers use big endian addressing,
# while Intel-based Macs use little-endian addressing. Most non-Intel-based computers are big endian.
# see: https://developer.apple.com/library/content/documentation/DeviceDrivers/Conceptual/WritingPCIDrivers/endianness/endianness.html


# we will need to specify endianness, otherwise:
# ..
# obspy.io.segy.segy.SEGYError: Unable to determine the endianness of the file. Please specify it.
# ..
endian = '<' # little endian '<' or big endian '>'

# Seismic Unix
if read_su:
    # reads as stream
    st = read(file_in,format='SU',byteorder=endian)

# SEG-Y
if read_segy:
    # reads as stream
    st = read(file_in,format='SEGY',byteorder=endian)

    ## more detailed call:
    #from obspy.io.segy.segy import _read_segy
    #import obspy.io as io

    # By default trace data are read into memory, but this may be impractical for very large datasets
    # see: https://docs.obspy.org/packages/obspy.io.segy.html
    #
    #segy = io.segy.segy._read_segy(file_in, headonly=True,endian=endian)
    #print(len(segy.traces[0].data))

    ## reading w/ data
    #segy = io.segy.segy._read_segy(file_in,endian=endian)
    #print(segy)

    ## reads as stream
    #st = io.segy.core._read_segy(file_in,byteorder=endian)

# output info
print("traces:")
print(st)
print("")

# timing
# how long it takes to determine minimum and maximum values in an array
from timeit import default_timer

def minmax_np(xy):
    # uses numpy internal functions
    minval = np.min(xy[:,1])
    maxval = np.max(xy[:,1])
    return (minval,maxval)

def minmax(xy):
    # combines loop over all elements to determine min/max at the same time
    minval = xy[0,1]
    maxval = xy[0,1]
    for f1 in xy[1:,1]:
        if f1 > maxval:
            maxval = f1
        elif f1 < minval:
            minval = f1
    return (minval, maxval)

def minmax2(xy):
    # combines loop over all elements and takes steps of 2 to minimize number of comparisons
    minval = xy[0,1]
    maxval = xy[0,1]
    if len(xy[:,1]) % 2 == 0:
      istart = 0
    else:
      istart = 1
    for ii in xrange(0, len(xy[:,1]), 2):
        f1 = xy[ii  ,1]
        f2 = xy[ii+1,1]
        if (f1 < f2):
            if f1 < minval: minval = f1
            if f2 > maxval: maxval = f2
        else:
            if f2 < minval: minval = f2
            if f1 > maxval: maxval = f1
    return (minval, maxval)


# direct looping without prior reading in of stream...
#for tr in obspy.io.segy.segy.iread_su('OUTPUT_FILES/0_dz_SU',endian='<'):

# prints traces
# see: https://docs.obspy.org/tutorial/code_snippets/export_seismograms_to_ascii.html

for i,tr in enumerate(st):
    # file output
    filename = str.format("%s_%d" % (out_file, i))
    print(filename)

    # length
    length = len(tr.data)
    t0 = 0.0
    dt = tr.stats.delta
    print("  length = ",length," dt = ",dt)

    # SPECFEM won't fill station and channel header infos
    # station name
    if len(tr.stats.station) > 0:
        sta = tr.stats.station
    else:
        sta = 'A' + str(i)                    # A0, A1, ..
    # channel name
    if len(tr.stats.channel) > 0:
        cha = tr.stats.channel
    else:
        # using file name for a guess
        basename = os.path.basename(file_in)  # 0_dz_SU
        comps = basename.upper().split("_")   # ['0','DZ','SU']
        if len(comps) > 1:
            if len(comps[1]) > 1:
                cha = comps[1][-1:]           # 'Z'
            else:
                cha = '*'
        else:
            cha = '*'

    # file output
    f = open(filename, "w")

    # file header
    f.write("# STATION %s\n" % (sta))
    f.write("# CHANNEL %s\n" % (cha))
    f.write("# START_TIME %s\n" % (str(tr.stats.starttime)))
    f.write("# SAMP_FREQ %f\n" % (tr.stats.sampling_rate))
    f.write("# NDAT %d\n" % (tr.stats.npts))

    # data section
    # fills numpy array
    xy = np.empty(shape=(0,2))
    for ii in range(length):
        time = ii*dt + t0
        xy = np.append(xy,[[time,tr.data[ii]]],axis=0)

    # data column
    #print(xy[:,1])
    #print(xy[0,1],xy[1,1],xy[2,1],xy[3,1])

    # saves as ascii
    np.savetxt(f, xy, fmt="%f\t%f")
    f.close()

    # statistics
    #start = default_timer()
    minval,maxval = minmax_np(xy)   # 7.414818e-05 s
    #minval,maxval = minmax(xy)     # 6.389618e-04 s
    #minval,maxval = minmax2(xy)    # 1.811981e-03 s
    #end = default_timer()
    #print("  min/max = %f / %f  (timing: %e s)" %(minval,maxval,end-start))
    print("  min/max = %f / %f" %(minval,maxval))
    print("  done")

# figure output
#st.plot()