Go to the documentation of this file.00001
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]))