latlon_to_utm.py
Go to the documentation of this file.
00001 #!/usr/bin/env python
00002 
00003 import pyproj
00004 import math
00005 import sys
00006 from optparse import OptionParser
00007 from math import atan2, cos, sin, sqrt
00008 
00009 
00010 def pythagoras(a, b):
00011     return sqrt(a**2 + b**2)
00012 
00013 
00014 class WGS84():
00015     A = 6378137.0
00016     recf = 298.257223563
00017     B = A * (recf - 1.0) / recf
00018     eSquared = ((A * A) - (B * B)) / (A * A)
00019     e2Squared = ((A * A) - (B * B)) / (B * B)
00020     tn = (A - B) / (A + B)
00021     ap = A * (1.0 - tn + 5.0 * ((tn * tn) - (tn * tn * tn)) / 4.0 + 81.0 * ((tn * tn * tn * tn) - (tn * tn * tn * tn * tn)) / 64.0)
00022     bp = 3.0 * A * (tn - (tn * tn) + 7.0 * ((tn * tn * tn) - (tn * tn * tn * tn)) / 8.0 + 55.0 * (tn * tn * tn * tn * tn) / 64.0) / 2.0
00023     cp = 15.0 * A * ((tn * tn) - (tn * tn * tn) + 3.0 * ((tn * tn * tn * tn) - (tn * tn * tn * tn * tn)) / 4.0) / 16.0
00024     dp = 35.0 * A * ((tn * tn * tn) - (tn * tn * tn * tn) + 11.0 * (tn * tn * tn * tn * tn) / 16.0) / 48.0
00025     ep = 315.0 * A * ((tn * tn * tn * tn) - (tn * tn * tn * tn * tn)) / 512.0
00026 
00027 
00028 def ecef_to_utm(x, y, z):
00029     xyz = (x, y, z)
00030     a2 = WGS84.A * WGS84.A
00031     b2 = WGS84.B * WGS84.B
00032     p = pythagoras(xyz[0], xyz[1])
00033     theta = atan2((xyz[2] * WGS84.A) , (p * WGS84.B))
00034     esq = 1.0 - (b2 / a2)
00035     epsq = (a2 / b2 - 1.0)
00036     sth = sin(theta)
00037     sth3 = sth**3
00038     cth = cos(theta)
00039     cth3 = cth**3
00040     lat = atan2((xyz[2] + epsq * WGS84.B * sth3) , (p - esq * WGS84.A * cth3))
00041     if (lat > (math.pi/2)):
00042         lat = math.pi - lat
00043     lon = atan2(xyz[1],xyz[0])
00044     hgt = (p / cos(lat)) - (a2 / sqrt(a2 * cos(lat) * cos(lat) +  b2 * sin(lat) * sin(lat)))
00045     lat = lat * (180/math.pi)
00046     lon = lon * (180/math.pi)
00047     return latlon_to_utm(lat, lon)
00048 
00049 
00050 def latlon_to_utm(latitude, longitude, height=0):
00051     """calculate UTM projection"""
00052     lon_6 = math.floor(longitude / 6.0)
00053     utmXZone = (longitude <= 0.0) and (30 + lon_6) or (31 + lon_6)
00054     p = pyproj.Proj(proj='utm', zone=utmXZone, ellps='WGS84')
00055     easting, northing = p(longitude, latitude)
00056     return (easting, northing)
00057 
00058 if __name__ == "__main__":
00059     parser = OptionParser()
00060     parser.add_option('--latitude', help='Latitude value in float degrees [-90.0:90]')
00061     parser.add_option('--longitude', help='Longitude value in float degrees [-180.0:180]')
00062     parser.add_option('-x', help='ECEF-X')
00063     parser.add_option('-y', help='ECEF-Y')
00064     parser.add_option('-z', help='ECEF-Z')
00065     (options, args) = parser.parse_args()
00066     if options.latitude and options.longitude:
00067         en = latlon_to_utm(float(options.latitude), float(options.longitude))
00068     elif options.x and options.y and options.z:
00069         en = ecef_to_utm(float(options.x), float(options.y), float(options.z))
00070     else:
00071         print("Select either latitude/longitude or x, y, z")
00072         sys.exit(1)
00073     print("easting northing: %f %f" % (en[0], en[1]))


rtklib
Author(s): Manuel Stahl
autogenerated on Mon Oct 6 2014 00:14:19