#! /usr/bin/env python
#--------------------------------------------------------------------------------------
# This script will convert all *.sgy files in a directory to Seismic Unix (SU) format,
# and from these files export the first four channels of every Shotpoint (SP) to
# ASCII files, for subsequent processing with e.g. Matlab or Octave.
#
# University of Bergen, Norway
# Dept. of Earth Science
# 19 April 2006
# O. Meyer
#--------------------------------------------------------------------------------------

from os import system,listdir,popen
from sys import *
from string import *
import time

# Define where SEGY, SU and ASCII files will live
data_dir  = '../data/'
SU_dir    =  data_dir + 'SU/'
ASCII_dir =  SU_dir   + 'ASCII/'

# Crate SU and ASCII directories
system('mkdir ' + SU_dir)
system('mkdir ' + ASCII_dir)

# Convert to SU format; strip 'sgy' file extension and add 'su'
print 'Converting all SEGY files to SU format ...'
system('for FILE in $(ls '+data_dir+'*.sgy); do FILE_A=${FILE##*/};\
        FILE_SU=${FILE_A%.*}.su; echo $FILE_SU; \
		segyread tape=$FILE conv=1 endian=0 | \
		segyclean > '+SU_dir+'$FILE_SU; done')
# 'FILE_A=${FILE##*/}' strips off directory part of file name, 'FILE_SU=${FILE_A%.*}.su' strips off '.sgy' extension and adds '.su'

# Now process each SU file: Find min/max SP in each file, export the four first channels to ASCII
filelist = listdir(SU_dir)
for file in filelist:
	if '.su' in SU_dir+file:
		print '-'*60
		print 'Now processing:', SU_dir+file 
		surange_list = popen('surange < ' + SU_dir+file)	# Use 'surange' command to get list of trace headers in use and their min/max values
		for line in surange_list.readlines():
			s = line.split()
			
			if 'fldr' in line:								# First get min/max fldr range from 'surange' command .......
				fldr_min = int(s[1])
				try:
				    fldr_max = int(s[2])
				except:
				    fldr_max = fldr_min
				print '    SP range: ['+str(fldr_min)+'..'+str(fldr_max)+']'
				
			if 'tracf' in line:								# Now get tracf min/max values
				tracf_min = int(s[1])
				try:
					tracf_max = int(s[2])
				except:
					tracf_max = tracf_min
					
				ch_extension = '-BPF-5-10-200-400-ch1+2+3+4.txt'			# Get 4 channels if possible, otherwise 3
				if tracf_max < 4:		
					tracf_max = 3
					ch_extension = '-BPF-5-10-200-400-ch1+2+3.txt'
				else:
					tracf_max = 4
				
			if 'ns' in line:								# Now get no of samples. Will use this later when we sort by traces
				ns = int(s[1])
			
				SP = fldr_min								# 'ns' is the last of these parametersin the 'surange' command; start processing ...
				# Loop over all Shotpoints (SP) in this file ...
				while SP <= fldr_max:
					filename = ASCII_dir+file+'.SP'+str(SP)+ch_extension		# Filename will be like '1166.su.SP_1166_ch_1+2+3+4.txt'
					print '    Writing '+filename
					# Here is the ASCII exporting; beware of the 'n1=4' option to the 'b2a' command if you make modifications later
					# The '2>' will get rid of status message from suwind that would mess up our listing
					system('rm temp.ascii')
					system('suwind key=fldr min='+str(SP)+' max='+str(SP)+' < '+SU_dir+file+' 2> /dev/null |\
					  suwind key=tracf min=1 max='+str(tracf_max)+' 2> /dev/null |\
					  sufilter f=5,10,200,400 | sustrip outpar=/dev/null | b2a n1=1 outpar=/dev/null > temp.ascii')
					  
					fd = open('./temp.ascii','r')		# This file holds 1 sample per line
					
					traces = []
					for k in range(tracf_max):
						trace = []
						for n in range(ns):
							t = fd.readline()
							trace.append(t)
						traces.append(trace)
					fd.close()
					
					fd = open(filename,'w')
					for n in range(ns):
						t = ''
						for k in range(tracf_max):
							s = traces[k][n]
							t = s[:-2]					# Lines in temp.ascii is terminated by SPACE + CR
							fd.write(t)
							#print t
							if (k <> (tracf_max-1)):
								fd.write(' ')
							else:
								fd.write('\n')

					fd.close()
						
					  
					SP += 1