#! /usr/bin/env python
# -*- coding: cp1252 -*-
'''
***************************************************************************************
                      S M O L A - 2 0 0 6   S U R V E Y
             
DESCRIPTION:
   This Python script tries to match SP and any one-minute position that is very close.
   Compare SP time and time of these matching one/minute positions.

   We must take into account that shotpoint log positions are given as UTM
   Easting/Northing Zone 31N, and interpolated data sets are given as UTM Zone 32N.
   So we must convert from UTM Zone 31 to Zone 32.

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.

VERSION:

    Ver  Date         By         Description
    ----------------------------------------------------------------------------------
    0.1  15 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, sleep
from cmath import *

#-----  Constants
UTM_Zone          = '31N'
WGS84             = 23                  # Datum constant, see LatLongUTMconversion file.
MAX_DISTANCE      = 1.0


#-----  Directory, file names
DataDirectory     = './logs'
ShotPointFile     = 'smoela_log.txt'
OneMinutePosFile  = 'toktlogger-24-25-okt-2006.txt'



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

#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#>>>>>>    1 . s t   s t a g e   <<<<<<<
#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

#---- Loop through all lines in EIVA shotpoint file
#---- Extract Shotpoint number, Easting/Northing; transform E/N into neighbouring UTM Zone

f = open(join(DataDirectory, ShotPointFile))

Eiva_Dict  = {}

for line in f:
    if line[0] == '0':
        s = line.split()        
        EventNo   = s[0]
        Date      = s[1]
        Time      = s[2]
        Easting   = s[3]
        Northing  = s[4]

        #--- Shotpoint Easting, Northing values in UTM Zone 31.
        #--- Transform to Zone 32, via geograpical coordinates.
        
        (lat, lon) = LatLongUTMconversion.UTMtoLL(WGS84, float(Northing), float(Easting), UTM_Zone)
        (z, e, n)  = LatLongUTMconversion.LLtoUTM(WGS84, lat, lon)
        #print EventNo, Easting, Northing, lat, lon, z, e, n

        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_Dict[EventNo]    = (e, n, lat, lon, EivaTimeStamp)

f.close()


#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#>>>>>>    2 . n d   s t a g e   <<<<<<<
#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

#---- Loop through all lines in one-minute position file. Extract Time and Position

f = open(join(DataDirectory, OneMinutePosFile))
OneMinTimePos   = []

for line in f:                           #<<<<<<-- Loop throgh all lines in input file
    s = line.split(',')
    Day       = s[2].split('.')[0]
    Month     = s[2].split('.')[1]
    Year      = s[2].split('.')[2]
    Hour      = s[3][0:2]
    Minute    = s[3][2:4]
    Seconds   = s[3][4:6]
    Microsec  = int(s[3].split('.')[1]) * 10000
    Latitude  = s[12].split()[0]    # s[12] holds string like '6324.99468 N'
    Longitude = s[13].split()[0]    # s[13] holds string like '00706.58098 E'

    #--- Convert lat/long to decimal degrees, and then to UTM Easting/Northing
    Latitude  = int(Latitude[0:2])  + (float(Latitude[2:])/60.0)
    Longitude = int(Longitude[0:3]) + (float(Longitude[3:])/60.0)
    (z, e, n) = LatLongUTMconversion.LLtoUTM(WGS84, Latitude, Longitude)
    
    #--- Create Datetime object from current time values
    #print '%s-%s-%s  %02s:%02s:%02s.%06d' % (Day, Month, Year, Hour, Minute, Seconds, Microsec)
    TimeStamp = datetime(   int(Year),
                            int(Month),
                            int(Day),
                            int(Hour),
                            int(Minute),
                            int(Seconds),
                            Microsec
                            )

    OneMinTimePos.append((TimeStamp, e, n))

f.close()



#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#>>>>>>    3 . r d   s t a g e   <<<<<<<
#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

SP_keys = Eiva_Dict.keys()
SP_keys.sort()

for SP in SP_keys:              #<<<--- Loop through every shotpoint
    (EivaEasting, EivaNorthing, EivaLat, EivaLong, EivaTimeStamp) = Eiva_Dict[SP]
    cShotpointPosition = EivaEasting + 1j*EivaNorthing

    for k in range(len(OneMinTimePos)):
        (TimeStamp, e, n) = OneMinTimePos[k]
        #print '::', k, TimeStamp, e, n
        cOneMinPosition = e + 1j*n
        cDeltaPos = cShotpointPosition - cOneMinPosition
        TimeDelta = EivaTimeStamp - TimeStamp
        if abs(cDeltaPos) < MAX_DISTANCE:
            print "%s %.2f %.2f" % (SP, EivaEasting, EivaNorthing),EivaTimeStamp, "%.2f %.2f" % (e, n),TimeStamp, "%.2f" % abs(cDeltaPos),TimeDelta

    
print 'Done.'