#!/usr/bin/env python # calculating ship speed GOS-2008 gravity data #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ # LAT/LONG to/from UTM library, from # http://www.pygps.org/#LatLongUTMconversion #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ import cmath, os # Lat Long - UTM, UTM - Lat Long conversions from math import pi, sin, cos, tan, sqrt #LatLong- UTM conversion..h #definitions for lat/long to UTM and UTM to lat/lng conversions #include _deg2rad = pi / 180.0 _rad2deg = 180.0 / pi _EquatorialRadius = 2 _eccentricitySquared = 3 _ellipsoid = [ # id, Ellipsoid name, Equatorial Radius, square of eccentricity # first once is a placeholder only, To allow array indices to match id numbers [ -1, "Placeholder", 0, 0], [ 1, "Airy", 6377563, 0.00667054], [ 2, "Australian National", 6378160, 0.006694542], [ 3, "Bessel 1841", 6377397, 0.006674372], [ 4, "Bessel 1841 (Nambia] ", 6377484, 0.006674372], [ 5, "Clarke 1866", 6378206, 0.006768658], [ 6, "Clarke 1880", 6378249, 0.006803511], [ 7, "Everest", 6377276, 0.006637847], [ 8, "Fischer 1960 (Mercury] ", 6378166, 0.006693422], [ 9, "Fischer 1968", 6378150, 0.006693422], [ 10, "GRS 1967", 6378160, 0.006694605], [ 11, "GRS 1980", 6378137, 0.00669438], [ 12, "Helmert 1906", 6378200, 0.006693422], [ 13, "Hough", 6378270, 0.00672267], [ 14, "International", 6378388, 0.00672267], [ 15, "Krassovsky", 6378245, 0.006693422], [ 16, "Modified Airy", 6377340, 0.00667054], [ 17, "Modified Everest", 6377304, 0.006637847], [ 18, "Modified Fischer 1960", 6378155, 0.006693422], [ 19, "South American 1969", 6378160, 0.006694542], [ 20, "WGS 60", 6378165, 0.006693422], [ 21, "WGS 66", 6378145, 0.006694542], [ 22, "WGS-72", 6378135, 0.006694318], [ 23, "WGS-84", 6378137, 0.00669438] ] #Reference ellipsoids derived from Peter H. Dana's website- #http://www.utexas.edu/depts/grg/gcraft/notes/datum/elist.html #Department of Geography, University of Texas at Austin #Internet: pdana@mail.utexas.edu #3/22/95 #Source #Defense Mapping Agency. 1987b. DMA Technical Report: Supplement to Department of Defense World Geodetic System #1984 Technical Report. Part I and II. Washington, DC: Defense Mapping Agency #def LLtoUTM(int ReferenceEllipsoid, const double Lat, const double Long, # double &UTMNorthing, double &UTMEasting, char* UTMZone) def LLtoUTM(ReferenceEllipsoid, Lat, Long, zone = None): '''converts lat/long to UTM coords. Equations from USGS Bulletin 1532 East Longitudes are positive, West longitudes are negative. North latitudes are positive, South latitudes are negative Lat and Long are in decimal degrees Written by Chuck Gantz- chuck.gantz at globalstar.com''' a = _ellipsoid[ReferenceEllipsoid][_EquatorialRadius] eccSquared = _ellipsoid[ReferenceEllipsoid][_eccentricitySquared] k0 = 0.9996 #Make sure the longitude is between -180.00 .. 179.9 LongTemp = (Long+180)-int((Long+180)/360)*360-180 # -180.00 .. 179.9 LatRad = Lat*_deg2rad LongRad = LongTemp*_deg2rad if zone is None: ZoneNumber = int((LongTemp + 180)/6) + 1 else: ZoneNumber = zone if Lat >= 56.0 and Lat < 64.0 and LongTemp >= 3.0 and LongTemp < 12.0: ZoneNumber = 32 # Special zones for Svalbard if Lat >= 72.0 and Lat < 84.0: if LongTemp >= 0.0 and LongTemp < 9.0:ZoneNumber = 31 elif LongTemp >= 9.0 and LongTemp < 21.0: ZoneNumber = 33 elif LongTemp >= 21.0 and LongTemp < 33.0: ZoneNumber = 35 elif LongTemp >= 33.0 and LongTemp < 42.0: ZoneNumber = 37 LongOrigin = (ZoneNumber - 1)*6 - 180 + 3 #+3 puts origin in middle of zone LongOriginRad = LongOrigin * _deg2rad #compute the UTM Zone from the latitude and longitude UTMZone = "%d%c" % (ZoneNumber, _UTMLetterDesignator(Lat)) eccPrimeSquared = (eccSquared)/(1-eccSquared) N = a/sqrt(1-eccSquared*sin(LatRad)*sin(LatRad)) T = tan(LatRad)*tan(LatRad) C = eccPrimeSquared*cos(LatRad)*cos(LatRad) A = cos(LatRad)*(LongRad-LongOriginRad) M = a*((1 - eccSquared/4 - 3*eccSquared*eccSquared/64 - 5*eccSquared*eccSquared*eccSquared/256)*LatRad - (3*eccSquared/8 + 3*eccSquared*eccSquared/32 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(2*LatRad) + (15*eccSquared*eccSquared/256 + 45*eccSquared*eccSquared*eccSquared/1024)*sin(4*LatRad) - (35*eccSquared*eccSquared*eccSquared/3072)*sin(6*LatRad)) UTMEasting = (k0*N*(A+(1-T+C)*A*A*A/6 + (5-18*T+T*T+72*C-58*eccPrimeSquared)*A*A*A*A*A/120) + 500000.0) UTMNorthing = (k0*(M+N*tan(LatRad)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24 + (61 -58*T +T*T +600*C -330*eccPrimeSquared)*A*A*A*A*A*A/720))) if Lat < 0: UTMNorthing = UTMNorthing + 10000000.0; #10000000 meter offset for southern hemisphere return (UTMZone, UTMEasting, UTMNorthing) def _UTMLetterDesignator(Lat): '''This routine determines the correct UTM letter designator for the given latitude returns 'Z' if latitude is outside the UTM limits of 84N to 80S Written by Chuck Gantz- chuck.gantz@globalstar.com''' if 84 >= Lat >= 72: return 'X' elif 72 > Lat >= 64: return 'W' elif 64 > Lat >= 56: return 'V' elif 56 > Lat >= 48: return 'U' elif 48 > Lat >= 40: return 'T' elif 40 > Lat >= 32: return 'S' elif 32 > Lat >= 24: return 'R' elif 24 > Lat >= 16: return 'Q' elif 16 > Lat >= 8: return 'P' elif 8 > Lat >= 0: return 'N' elif 0 > Lat >= -8: return 'M' elif -8> Lat >= -16: return 'L' elif -16 > Lat >= -24: return 'K' elif -24 > Lat >= -32: return 'J' elif -32 > Lat >= -40: return 'H' elif -40 > Lat >= -48: return 'G' elif -48 > Lat >= -56: return 'F' elif -56 > Lat >= -64: return 'E' elif -64 > Lat >= -72: return 'D' elif -72 > Lat >= -80: return 'C' else: return 'Z' # if the Latitude is outside the UTM limits #void UTMtoLL(int ReferenceEllipsoid, const double UTMNorthing, const double UTMEasting, const char* UTMZone, # double& Lat, double& Long ) def UTMtoLL(ReferenceEllipsoid, northing, easting, zone): """converts UTM coords to lat/long. Equations from USGS Bulletin 1532 East Longitudes are positive, West longitudes are negative. North latitudes are positive, South latitudes are negative Lat and Long are in decimal degrees. Written by Chuck Gantz- chuck.gantz@globalstar.com Converted to Python by Russ Nelson """ k0 = 0.9996 a = _ellipsoid[ReferenceEllipsoid][_EquatorialRadius] eccSquared = _ellipsoid[ReferenceEllipsoid][_eccentricitySquared] e1 = (1-sqrt(1-eccSquared))/(1+sqrt(1-eccSquared)) #NorthernHemisphere; //1 for northern hemispher, 0 for southern x = easting - 500000.0 #remove 500,000 meter offset for longitude y = northing ZoneLetter = zone[-1] ZoneNumber = int(zone[:-1]) if ZoneLetter >= 'N': NorthernHemisphere = 1 # point is in northern hemisphere else: NorthernHemisphere = 0 # point is in southern hemisphere y -= 10000000.0 # remove 10,000,000 meter offset used for southern hemisphere LongOrigin = (ZoneNumber - 1)*6 - 180 + 3 # +3 puts origin in middle of zone eccPrimeSquared = (eccSquared)/(1-eccSquared) M = y / k0 mu = M/(a*(1-eccSquared/4-3*eccSquared*eccSquared/64-5*eccSquared*eccSquared*eccSquared/256)) phi1Rad = (mu + (3*e1/2-27*e1*e1*e1/32)*sin(2*mu) + (21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(4*mu) +(151*e1*e1*e1/96)*sin(6*mu)) phi1 = phi1Rad*_rad2deg; N1 = a/sqrt(1-eccSquared*sin(phi1Rad)*sin(phi1Rad)) T1 = tan(phi1Rad)*tan(phi1Rad) C1 = eccPrimeSquared*cos(phi1Rad)*cos(phi1Rad) R1 = a*(1-eccSquared)/pow(1-eccSquared*sin(phi1Rad)*sin(phi1Rad), 1.5) D = x/(N1*k0) Lat = phi1Rad - (N1*tan(phi1Rad)/R1)*(D*D/2-(5+3*T1+10*C1-4*C1*C1-9*eccPrimeSquared)*D*D*D*D/24 +(61+90*T1+298*C1+45*T1*T1-252*eccPrimeSquared-3*C1*C1)*D*D*D*D*D*D/720) Lat = Lat * _rad2deg Long = (D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*eccPrimeSquared+24*T1*T1) *D*D*D*D*D/120)/cos(phi1Rad) Long = LongOrigin + Long * _rad2deg return (Lat, Long) #----------------------------------------------------------------------------------------------------------------------- # S T A R T O F P R O G R A M #----------------------------------------------------------------------------------------------------------------------- ''' Lines in input file can look like this: 1215916740,13.07.2008,023900,73.594672426,7.692507784,2500.9,10079.61 1215916750,13.07.2008,023910,73.594308518,7.691386309,2502.0,10079.45 ''' # 1knot = 1852 m per hour => 1852/3600 = 0.514444444 m/s KnotFactor = 1852.0/3600.0 FILE = "sars-2008-possible-lines.txt" f = open(FILE) StartNewSegment = True NewSegment = '' PresentSegment = '' FirstPositionInSegment = False PreviousTimestamp = '' #print "#Timestamp[Unix], Lat, Long, UTMzone, Easting, Northing, Distance, dTime, Speed[m/s], Speed[knots]" for line in f: if len(line) <> 0: if "LINE" in line: try: g.close() h = os.popen('gnuplot' ,'w') print >>h, "set terminal png large size 1200,1024; set out '%s'" % PNGFileName print >>h, "set xdata time" print >>h, "set yrange [0:12]" print >>h, "set timefmt '%s'" print >>h, "set format x '%j_%H:%M'" print >>h, "set title 'Gravity data, GOS-2008. Plot of ship speed for %s'" % NewSegment print >>h, "set grid; set xlabel 'Time [UTC, day of year and time]'; set ylabel 'Ship Speed [kt]'" print >>h, "plot '%s' using 1:10 title 'Ship speed' with lines lw 4" % TextFileName #print >>h, "replot" h.flush() h.close() except: pass NewSegment = line.strip() TextFileName = "gravity-GOS2008-line-" + ("%02d" % int(NewSegment.split('LINE')[1])) PNGFileName = TextFileName + ".png" TextFileName = TextFileName + ".txt" g = open(TextFileName, "w") g.write("#Timestamp[Unix], Lat, Long, UTMzone, Easting, Northing, Distance, dTime, Speed[m/s], Speed[knots]\n") if PresentSegment <> NewSegment: print "Processing:", NewSegment print "--- Text output file:", TextFileName print "--- Plot image file:", PNGFileName PresentSegment = NewSegment FirstPositionInSegment = True if line[0].isdigit(): #print line.strip() s = line.split(',') Timestamp = s[0] Lat = s[3] Long = s[4] (z, e, n) = LLtoUTM(23, float(Lat), float(Long)) if FirstPositionInSegment: cPreviousPosition = e + 1j*n PreviousTimestamp = Timestamp FirstPositionInSegment = False print else: cThisPosition = e + 1j*n cDeltaPosition = cThisPosition - cPreviousPosition Distance = abs(cDeltaPosition) dTime = int(Timestamp) - int(PreviousTimestamp) cPreviousPosition = cThisPosition PreviousTimestamp = Timestamp Speed = Distance/dTime SpeedKt = Speed/KnotFactor OutputText = "%s %s %s %s %1.1f %1.1f %1.1f %i %1.2f %1.2f\n" % (Timestamp, Lat, Long, z, e, n, Distance, dTime, Speed, SpeedKt) #print Lat, Long, #print "%s %1.1f %1.1f" % (z, e, n), #print "%1.1f %i %1.2f %1.2f" % (Distance, dTime, Speed, SpeedKt) #print OutputText, g.write(OutputText) f.close() g.close() #--- For last line segment - should be improved ..... h = os.popen('gnuplot' ,'w') print >>h, "set terminal png large size 1200,1024; set out '%s'" % PNGFileName print >>h, "set xdata time" print >>h, "set yrange [0:12]" print >>h, "set timefmt '%s'" print >>h, "set format x '%j_%H:%M'" print >>h, "set title 'Gravity data, GOS-2008. Plot of ship speed for %s'" % NewSegment print >>h, "set grid; set xlabel 'Time [UTC, day of year and time]'; set ylabel 'Ship Speed [kt]'" print >>h, "plot '%s' using 1:10 title 'Ship speed' with lines lw 4" % TextFileName #print >>h, "replot" h.flush() h.close()