#! /usr/bin/env python
# -*- coding: cp1252 -*-
'''
***************************************************************************************
                      S V A L E X - 2 0 0 7   S U R V E Y

                    U K O O A   F I L E   G E N E R A T I O N
             

DESCRIPTION:
   Script applies an offset to easting, northing to adjust for distance between
   seismic source and position reference point. 

   Then UKOOA and plot files are generated.

INPUT PARAMETERS:
   No. 1: Input file
          Data column description: See extraction in first section.
   No. 2: UTM Zone
   No. 3: Line name
   No. 4: Offset distance
   

OUTPUT FILES:
    Will be placed in same directory as input file.

UKOOA HEADER, UTM ZONE, OFFSET DISTANCE
    See 'Constant' section below.
    The UKOOA automatically includes values from other constants.

EXTERNAL LIBRARY USED:
    UTM to Lat/Long conversion - download library from:
    http://www.pygps.org/#LatLongUTMconversion
    Either:
    a) On Linux: Unpack it, and run "python setup.py install" as root.
       In order to build Python modules you must first install 'python-devel'.
    b) Windows/Linux: put "LatLongUTMconversion.py" in the same directory as
       this file.

PYTHON ON WINDOWS:
    Script was developed on Windows version of Python. It should run on
    Unix/Linux without any changes. Download Python for Windows from
    http://www.python.org
    N O T E :  In UKOOA file, line termination will be CR+LF on Win, and LF
    on Linux/Unix.
    
VERSION:

    Ver  Date         By         Description
    ----------------------------------------------------------------------------------
    0.1  29 Aug 2007  O. Meyer   Initial version
    ----------------------------------------------------------------------------------

                               Dept. of Earth Science
                                University of Bergen
                                    N O R W A Y

***************************************************************************************
'''

import LatLongUTMconversion
from   string   import split
from   os.path  import join
from   datetime import datetime, timedelta
from   time     import strptime, strftime, gmtime
from   cmath    import *
from   sys import *


#---- First check number of command line parameters
if len(argv[0:]) <> 5:		# One more that actual number, as name of script is in argv[0]
    print"Error: Need arguments: Navdata file, UTM Zone, Line Name, Source"
    exit(1)

#---- Get command line parameters:
#  No. 1: Input file
#         Data column description: See extraction in first section.
#   No. 2: UTM Zone
#   No. 3: Line name
#   No. 4: Offset distance

DATA_FILE     = argv[1]
UTM_Zone      = argv[2]
LINE_NAME     = argv[3]
SOURCE_OFFSET = argv[4]

#-----  Directory, file names

DATA_DIRECTORY    = './VAN_MIJEN_FJORD/EIVA'
#DATA_FILE         = 'SvalexVMF-07-01'
UKOOA_FILE        = DATA_FILE + '.UKOOA.txt'
PLOT_FILE         = DATA_FILE + '.plot.txt'

#-----  Constants

#LINE_NAME             = 'VM-01'         # Max 5 characters
#UTM_Zone              = '33X'
WGS84                 = 23              # Datum constant, see LatLongUTMconversion file.

#SOURCE_OFFSET         = 57.0            # Offset shotpoint positions by this amount.

FIRST_SP              = 1
LAST_SP               = 1086
CR_LF                 = '\x0d\x0a'

UKOOA_Header          = \
'''H0100 SURVEY AREA VAN MIJEN FJORD, SVALBARD, NORWAY
H0101 SURVEY DETAILS SVALEX-2007
H0102 VESSEL DETAILS R.V. HAAKON MOSBY
H0200 SURVEY DATE AUGUST 2007
H0201 FILE CREATED %s
H0202 FILE VERSION UKOOA P1/1990
H0300 CLIENT SVALEX-2007
H0400 GEOPHYSICAL CONTRACTOR UNIVERSITY OF BERGEN, NORWAY
H0500 POSITIONING CONTRACTOR N/A
H0600 POSITIONING PROCESSING O.M. UNIVERSITY OF BERGEN, DEPT OF EARTH SCIENCE
H0700 POSITIONING SYSTEM
H0700 SURVEY COMPUTER SOFTWARE NAVIPAC, EIVA, DENMARK
H0700 GPS SEAPATH 200 (PRIM), TRIMBLE 200D (SEC)
H0800 SHOTPOINT POSITION CENTRE OF SOURCE ARRAY
H1000 CLOCK TIME UTC
H0900 OFFSET SHIP SYSTEM TO SP %s M BEHIND
H1200 SPHEROID AS SURVEYED A,1/F WGS-84 Spheroid 6378137.000 298.2572236
H1300 SPHEROID AS PLOTTED A,1/F WGS-84 Spheroid 6378137.000 298.2572236
H1400 GEODETIC DATUM AS SURVEYED WGS-84
H1500 GEODETIC DATUM AS PLOTTED WGS-84
H1700 VERTICAL DATUM SL : ECHO SOUNDER
H1800 PROJECTION 001 UTM NORTHERN HEMISPHERE
H1900 ZONE %s NORTHERN HEMISPHERE
H2000 GRID UNIT 1 METRE
H2001 HEIGHT UNIT 1 METRE
H2002 ANGULAR UNITS 1 DEGREES
H2600 SOURCE 6 AIRGUNS, PRESSURE 138 BAR
H2600
H2600 00011111111112222222222333333333344444444445555555555666666666677777777778
H2600 78901234567890123456789012345678901234567890123456789012345678901234567890
H2600              <-SP-><---LAT--><--LONG---><--EAST-><-NORTH-><-DEP><D><TIME>-\n''' % (strftime("%d %b %Y %H:%M:%S UTC", gmtime()),
                                                                                         SOURCE_OFFSET,
                                                                                         UTM_Zone)



#-----------------------------------------------------------------------------------
#  -----     S U B R O U T I N E S
#-----------------------------------------------------------------------------------
def Degrees2DMS(desideg):            # Convert from desimal degrees to dddmmss.ss.
#-----------------------------------------------------------------------------------
    long_correction = float(desideg)			# These 4 lines will convert LONG in 0-360 range to 0-180 range
    if long_correction > 180.0:				# !!!! NOTE !!!! Take care of E/W 
        long_correction = 360.0 - long_correction
        desideg = str(long_correction)
    s = desideg.split('.')
    degree_main = s[0]
    degree_fraction = float(desideg) - float(degree_main)
    minute = str(degree_fraction * 60)
    minute_main = minute.split('.')[0]
    minute_fraction = '0.' + minute.split('.')[1]
    sec = float(minute_fraction) * 60
    second = "%5.2f" % sec

    return("%s%2d%s" % (degree_main, int(minute_main), second))



#---------------------------------------------------------------------
#-----            S t a r t   o f   p r o g r a m
#---------------------------------------------------------------------

#---- First check number of command line parameters
if len(argv[0:]) <> 6:		# One more that actual number, as name of script is in argv[0]
    print"Error: Need arguments: Navdata file, GPS timestamp file, nominal geometry file, UTM Zone, Line name"
    exit(1)

#---- Get command line parameters:
#  No. 1: Input file
#         Data column description: See extraction in first section.
#   No. 2: UTM Zone
#   No. 3: Line name
#   No. 4: Offset distance

DATA_FILE     = argv[1]
UTM_Zone      = argv[2]
LineName      = argv[3]
SOURCE_OFFSET = argv[4]




#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#>>>>>>       1 . s t   s e c t i o n       <<<<<<<
#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

#---- Loop through all lines in EIVA shotpoint file, and extract information.
#---- Convert UTM to Lat, Lon. Correct timestamp. Store in dictionary for later use.

f = open(join(DATA_DIRECTORY,DATA_FILE))

Eiva_Dictionary  = {}       # Will store all info extracted here. SP number used as key.
First_SP = 0

for line in f:
    if line[0] == '0':      # Assume SP is first, with leading zeros. Gets rid of comments, header this way.
        s = line.split()        
        Shotpoint = int(s[0])
        if First_SP == 0:
            First_SP = Shotpoint
        Date      = s[1]
        Time      = s[2]
        Easting   = float(s[3])
        Northing  = float(s[4])
        Depth     = float(s[12])

        #--- Convert e,n to lat,lon
        (Lat, Lon) = LatLongUTMconversion.UTMtoLL(WGS84, Easting, Northing, UTM_Zone)

        #---- Use date & time info to declare Python time object.
        #---- Using these time objects makes it easy to perform time adjustments.

        Year     = Date.split('.')[0]
        Month    = Date.split('.')[1]
        Day      = Date.split('.')[2]
        Hour     = Time.split(':')[0]
        Minute   = Time.split(':')[1]
        Second   = Time.split(':')[2].split('.')[0]
        Microsec = Time.split(':')[2].split('.')[1]

        EivaTimeStamp = datetime(int(Year),
                                 int(Month),
                                 int(Day),
                                 int(Hour),
                                 int(Minute),
                                 int(Second),
                                 (int(Microsec) * 1000)
                                 )

        #--- Store info in dictionary (with SP as key):    
        Eiva_Dictionary[Shotpoint]  = (Easting, Northing, EivaTimeStamp, Depth, Lat, Lon)
        #print SP_number, Date, Time, EivaTimeStamp, Easting, Northing, Lat, Lon, Depth      # Uncomment if you like watching scrolling text ...

f.close()



#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#>>>>>>           2 . n d   s e c t i o n         <<<<<<<
#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

#--  Offset adjustment to seismic source.
#--  Using vector calculus, let complex numbers represent vectors
#--  Offset adjustment it based on vectors: The vector between this
#--  and the previous SP is calculated. This vector is divided by
#--  its own length to obtain a unit length vector, and this in turn
#--  is multiplied by the offset distance. The resulting correction
#--  vector is subtracted from this shotpoint. First SP is treated
#--  differently, as there is no previous SP in that case; instead,
#--  use next SP to obtain correction vector.

f = open(join(DATA_DIRECTORY,PLOT_FILE), "w")
#-- Print header line in file
f.write("#SP, oldEasting, oldNorthing, oldLat, oldLon, newEasting, newNorthing, newLat, newLon\n")

SP_keys = Eiva_Dictionary.keys()
SP_keys.sort()

PositionDictionary = {}

for SP in SP_keys:              #<<<--- Loop through every shotpoint
    (ThisEasting, ThisNorthing, ThisTimeStamp, ThisDepth, ThisLat, ThisLon) = Eiva_Dictionary[SP]       # Retrieve info from dictionary, SP = key
    cThisPosition            = ThisEasting + 1j*ThisNorthing
    
    if SP != First_SP:
        (PrevEasting, PrevNorthing, PrevTimeStamp, PrevDepth, PrevLat, PrevLon) = Eiva_Dictionary[SP-1] 
        cPrevPosition        = PrevEasting + 1j*PrevNorthing
        cDeltaPosition       = cThisPosition - cPrevPosition
    else:   # This was the first SP, must treat it differently
        (NextEasting, NextNorthing, NextTimeStamp, NextDepth, NextLat, NextLon) = Eiva_Dictionary[SP+1]
        cNextPosition        = NextEasting + 1j*NextNorthing
        cDeltaPosition       = cNextPosition - cThisPosition
        
    cOffsetVector            = SOURCE_OFFSET * cDeltaPosition/abs(cDeltaPosition)
    cOffsetCorrectedPosition = cThisPosition - cOffsetVector

    #--- Transform new position to geograpical coordinates.
    e = cOffsetCorrectedPosition.real
    n = cOffsetCorrectedPosition.imag
    (NewLat, NewLon) = LatLongUTMconversion.UTMtoLL(WGS84, e, n, UTM_Zone)
    
    PositionDictionary[SP]   = (cOffsetCorrectedPosition, NewLat, NewLon)
    #print SP, cThisPosition, cOffsetVector, cOffsetCorrectedPosition       # Uncomment to see raw complex numbers in action
    # Plot original and offset adjusted positions
    f.write("%s %.2f %.2f %.10f %.10f %.2f %.2f %.10f %.10f\n" % (SP,
                                                                  cThisPosition.real,
                                                                  cThisPosition.imag,
                                                                  ThisLat,
                                                                  ThisLon,
                                                                  cOffsetCorrectedPosition.real,
                                                                  cOffsetCorrectedPosition.imag,
                                                                  NewLat,
                                                                  NewLon))
    
f.close()
print 'Plot file  : ' + PLOT_FILE


#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#>>>>>>           3 .  s e c t i o n              <<<<<<<
#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

# Write UKOOA P1/90 lines to file
# WARNING: Assumes Latitudes in Northern hemisphere, and Longitude in Eastern!

f = open(join(DATA_DIRECTORY,UKOOA_FILE), "w")

Buffer = UKOOA_Header

for SP in SP_keys:
    (OldEasting, OldNorthing, TimeStamp, Depth, OldLat, OldLon) = Eiva_Dictionary[SP]
    (cOffsetCorrectedPosition, Lat, Lon)                        = PositionDictionary[SP]
    Easting  = cOffsetCorrectedPosition.real
    Northing = cOffsetCorrectedPosition.imag
    
    time   = "%02d%02d%02d.%03.0f000" % (  TimeStamp.hour,
                                           TimeStamp.minute,
                                           TimeStamp.second,
                                           round(TimeStamp.microsecond/1000.0))

    ukooa  = 'S'                                # First char of record = S, meaning "Centre of source"
    if len(LINE_NAME) > 5:
        print 'ERROR: Line name:', LINE_NAME, 'too long - max 5 characters, correct ...'
        exit (1)
    ukooa += "%-5s"  % LINE_NAME	        # -5s means string, width=5, left justify
    ukooa += ' ' * 13
    ukooa += "%6d"   % SP		        # Event = SP; width = 6, right justified
    ukooa += "%9s"   % Degrees2DMS(str(Lat))    # Latitude
    ukooa += 'N'
    ukooa += "%10s"  % Degrees2DMS(str(Lon))    # Longitude
    ukooa += 'E'
    ukooa += "%9.1f" % Easting                  # Easting, 9 characters in total, one decimal point
    ukooa += "%9.1f" % Northing                 # Northing, --"--
    ukooa += "%6.1f" % Depth                    # Depth
    ukooa += "%3d"   % TimeStamp.timetuple()[7] # Day of year
    ukooa += time.split('.')[0]                 # Time
    ukooa += ' '                                # Spare field, width = 1
    ukooa += '\n'                               # This will generate CR+LF on Windows ....
              
    Buffer += ukooa
    if len(ukooa) > 82:
        print 'ERROR: UKOOA data line exceeding 82 characters'
        exit(1)
    #print ukooa,

Buffer += 'EOF\n'

f.write(Buffer)

f.close()
print 'UKOOA file : ' + UKOOA_FILE