#! /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 combines information from three files: - Shotpoint file - Interpolated Time/Easting, and Time/Northing The aim is to tie the position of each shotpoint to the closest position of the interpolated data sets, and read the time from the (interpolated) point. 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 14 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. GUNCO_TIME_OFFSET = 0.128 # In seconds. Add to EIVA shooting time, to account for GUNCO delay. BASE_TIME = datetime(2006,10,24,0,0,0,0) #----- Directory, file names DataDirectory = './logs' ShotPointFile = 'smoela_log.txt' TimeEastingFile = 'spline-time_easting.txt' TimeNorthingFile = 'spline-time_northing.txt' #--------------------------------------------------------------------- #----- S t a r t o f p r o g r a m #--------------------------------------------------------------------- #>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> #>>>>>> 1 . s t s t a g e <<<<<<< #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< #---- Loop throgh all lines in shotpoint file #---- Extract Shotpoint number, Easting/Northing; transform E/N into neighbouring UTM Zone f = open(join(DataDirectory, ShotPointFile)) SP_Pos = {} SP_PosGeo = {} SP_TimeDict = {} 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 #--- Store info in three dictionaries (all with SP as key): SP_Pos[EventNo] = (e, n) SP_PosGeo[EventNo] = (lat,lon) SP_TimeDict[EventNo] = (Date, Time) f.close() #>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> #>>>>>> 2 . n d s t a g e <<<<<<< #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< #---- Loop through all lines in Time/Easting and Time/Northing interpolated data files. #---- Combine information so we have a dictionary with Time as key, and (Easting, Northing) tuple #---- as value. fE = open(TimeEastingFile) fN = open(TimeNorthingFile) Time_Pos = {} for lineE in fE: s = lineE.split() TimeEasting = s[0] Easting = s[1] lineN = fN.readline() t = lineN.split() TimeNorthing = t[0] Northing = t[1] if TimeEasting != TimeNorthing: print 'Time mismatch in Time/Easting and Time/Northing files!' exit Time_Pos[TimeEasting] = (Easting, Northing) fE.close() fN.close() #>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> #>>>>>> 3 . r d s t a g e <<<<<<< #<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< #---- For each shotpoint position, find nearest position among interpolated #---- Time/Position values. #---- Uses complex math so that real and imaginary values correspond to Easting #---- and Northing, respectively. This facilitates distance calculations. #---- The search method is awful and also time consuming (for each SP check about #---- 130 000 or 1 300 000 intertpolated values - depending on which dataset is used, #---- the one with 100x interpolation, or 1000x interpolation). Will improve it later, #---- if there is time. #---- NOTE: The Survey logger was switched off in the first part of the line. #---- We have one-minute time/position data from SP=144, so just skip #---- earlier SP. Must use another correction method for these. Time_keys = Time_Pos.keys() Time_keys.sort() SP_keys = SP_Pos.keys() SP_keys.sort() for SP in SP_keys: #<<<--- Loop through every shotpoint if int(SP) > 143: # No interpolated values for SP less then 143, skip them. (SP_Easting, SP_Northing) = SP_Pos[SP] cShotpointPosition = SP_Easting + 1j*SP_Northing Dist = {} for T in Time_keys: #<<<--- For each shotpoint, loop through every interpolated position. (e, n) = Time_Pos[T] #<<<--- This method is awful - but it didn't take any programming time. cInterpolatedPosition = float(e) + 1j*float(n) #-- Calculate distance vector between SP and this interpolated point cDistance = cShotpointPosition - cInterpolatedPosition #-- Store length of distance vector in dictionary. Dist[abs(cDistance)] = T # Complex math is doing all the dirty work for us ... #-- Now we are done checking all interpolated points. Get point with minimum distance to SP location. k = min(Dist.keys()) # Keys are the different distances found. Value is Time. #-- So far Time has been relative to basetime at midnight on 24. October 2006. #-- We want 'true' time so add relative time to basetime. #-- Also add offset introduced by GUNCO (a constant) DeltaT = timedelta(0, float(Dist[k])) SP_Time = BASE_TIME + DeltaT + timedelta(0, GUNCO_TIME_OFFSET) Microsec = "%1.3f" % (SP_Time.microsecond/1000000.0) Microsec = Microsec.split('.')[1] (EivaDate, EivaTime) = SP_TimeDict[SP] EivaYear = EivaDate.split('.')[0] EivaMonth = EivaDate.split('.')[1] EivaDay = EivaDate.split('.')[2] EivaHour = EivaTime.split(':')[0] EivaMinute = EivaTime.split(':')[1] EivaSecond = EivaTime.split(':')[2].split('.')[0] EivaMicrosec = EivaTime.split(':')[2].split('.')[1] EivaTimeStamp = datetime(int(EivaYear), int(EivaMonth), int(EivaDay), int(EivaHour), int(EivaMinute), int(EivaSecond), (int(EivaMicrosec) * 1000) ) TimeDeviation = EivaTimeStamp - SP_Time #--------------------------------------------------------------------- #----- P r i n t o u t s e c t i o n #--------------------------------------------------------------------- # Parameter description: # # 1) SP number # 2) SP Date # 3) SP Time # 4) SP Easting # 5) SP Northing # 6) SP Lat # 7) SP Long # 8) Seconds since midnight 24 Oct 2006 # 9) Distance between SP and nearest interpolated point # 10) Easting of nearest interpolated point # 11) Northing of nearest interpolated point # 12) Difference between original time stamp (Eiva # computer local clock), and that of nearest interpolated point # #---------------------------------------------------------------------- print SP, print "%d.%02d.%02d %02d:%02d:%02d.%s" % (SP_Time.year, SP_Time.month, SP_Time.day, SP_Time.hour, SP_Time.minute, SP_Time.second, Microsec), print "%1.1f %1.1f" % (SP_Easting, SP_Northing), (lat,lon) = SP_PosGeo[SP] print "%1.10f %1.10f" % (lat, lon), print "%1.3f %1.1f" % (float(Dist[k]), k), # Time of nearest interpolated point, distance beween SP and this point. (e, n) = Time_Pos[Dist[k]] # Get position of this particular time value print "%1.1f %1.1f" % (float(e), float(n)), # Prints Easting, Northing of closest interpolated point. print "%1.2f" % (TimeDeviation.seconds + (TimeDeviation.microseconds/1000000.0)) print 'Done.'