Created
April 29, 2016 11:49
-
-
Save Alliages/33d50c5f3a75da0cb065698a5a011e52 to your computer and use it in GitHub Desktop.
Convert a (latitude, longitude) tuple into an UTM coordinate
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
""" | |
Authors | |
Bart van Andel <bavanandel@gmail.com> | |
Tobias Bieniek <Tobias.Bieniek@gmx.de> | |
Torstein I. Bø | |
License | |
Copyright (C) 2012 Tobias Bieniek <Tobias.Bieniek@gmx.de> | |
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so | |
https://github.com/Turbo87/utm | |
------ | |
usage : | |
import utm | |
location = utm.from_latlon(lat,lon) | |
print location[0] | |
print location[1] | |
------- | |
""" | |
import math | |
#from utm.error import OutOfRangeError | |
class OutOfRangeError(ValueError): | |
pass | |
__all__ = ['to_latlon', 'from_latlon'] | |
K0 = 0.9996 | |
E = 0.00669438 | |
E2 = E * E | |
E3 = E2 * E | |
E_P2 = E / (1.0 - E) | |
SQRT_E = math.sqrt(1 - E) | |
_E = (1 - SQRT_E) / (1 + SQRT_E) | |
_E2 = _E * _E | |
_E3 = _E2 * _E | |
_E4 = _E3 * _E | |
_E5 = _E4 * _E | |
M1 = (1 - E / 4 - 3 * E2 / 64 - 5 * E3 / 256) | |
M2 = (3 * E / 8 + 3 * E2 / 32 + 45 * E3 / 1024) | |
M3 = (15 * E2 / 256 + 45 * E3 / 1024) | |
M4 = (35 * E3 / 3072) | |
P2 = (3. / 2 * _E - 27. / 32 * _E3 + 269. / 512 * _E5) | |
P3 = (21. / 16 * _E2 - 55. / 32 * _E4) | |
P4 = (151. / 96 * _E3 - 417. / 128 * _E5) | |
P5 = (1097. / 512 * _E4) | |
R = 6378137 | |
ZONE_LETTERS = [ | |
(84, None), (72, 'X'), (64, 'W'), (56, 'V'), (48, 'U'), (40, 'T'), | |
(32, 'S'), (24, 'R'), (16, 'Q'), (8, 'P'), (0, 'N'), (-8, 'M'), (-16, 'L'), | |
(-24, 'K'), (-32, 'J'), (-40, 'H'), (-48, 'G'), (-56, 'F'), (-64, 'E'), | |
(-72, 'D'), (-80, 'C') | |
] | |
def to_latlon(easting, northing, zone_number, zone_letter=None, northern=None): | |
if not zone_letter and northern is None: | |
raise ValueError('either zone_letter or northern needs to be set') | |
elif zone_letter and northern is not None: | |
raise ValueError('set either zone_letter or northern, but not both') | |
if not 100000 <= easting < 1000000: | |
raise OutOfRangeError('easting out of range (must be between 100.000 m and 999.999 m)') | |
if not 0 <= northing <= 10000000: | |
raise OutOfRangeError('northing out of range (must be between 0 m and 10.000.000 m)') | |
if not 1 <= zone_number <= 60: | |
raise OutOfRangeError('zone number out of range (must be between 1 and 60)') | |
if zone_letter: | |
zone_letter = zone_letter.upper() | |
if not 'C' <= zone_letter <= 'X' or zone_letter in ['I', 'O']: | |
raise OutOfRangeError('zone letter out of range (must be between C and X)') | |
northern = (zone_letter >= 'N') | |
x = easting - 500000 | |
y = northing | |
if not northern: | |
y -= 10000000 | |
m = y / K0 | |
mu = m / (R * M1) | |
p_rad = (mu + | |
P2 * math.sin(2 * mu) + | |
P3 * math.sin(4 * mu) + | |
P4 * math.sin(6 * mu) + | |
P5 * math.sin(8 * mu)) | |
p_sin = math.sin(p_rad) | |
p_sin2 = p_sin * p_sin | |
p_cos = math.cos(p_rad) | |
p_tan = p_sin / p_cos | |
p_tan2 = p_tan * p_tan | |
p_tan4 = p_tan2 * p_tan2 | |
ep_sin = 1 - E * p_sin2 | |
ep_sin_sqrt = math.sqrt(1 - E * p_sin2) | |
n = R / ep_sin_sqrt | |
r = (1 - E) / ep_sin | |
c = _E * p_cos**2 | |
c2 = c * c | |
d = x / (n * K0) | |
d2 = d * d | |
d3 = d2 * d | |
d4 = d3 * d | |
d5 = d4 * d | |
d6 = d5 * d | |
latitude = (p_rad - (p_tan / r) * | |
(d2 / 2 - | |
d4 / 24 * (5 + 3 * p_tan2 + 10 * c - 4 * c2 - 9 * E_P2)) + | |
d6 / 720 * (61 + 90 * p_tan2 + 298 * c + 45 * p_tan4 - 252 * E_P2 - 3 * c2)) | |
longitude = (d - | |
d3 / 6 * (1 + 2 * p_tan2 + c) + | |
d5 / 120 * (5 - 2 * c + 28 * p_tan2 - 3 * c2 + 8 * E_P2 + 24 * p_tan4)) / p_cos | |
return (math.degrees(latitude), | |
math.degrees(longitude) + zone_number_to_central_longitude(zone_number)) | |
def from_latlon(latitude, longitude, force_zone_number=None): | |
if not -80.0 <= latitude <= 84.0: | |
raise OutOfRangeError('latitude out of range (must be between 80 deg S and 84 deg N)') | |
if not -180.0 <= longitude <= 180.0: | |
raise OutOfRangeError('longitude out of range (must be between 180 deg W and 180 deg E)') | |
lat_rad = math.radians(latitude) | |
lat_sin = math.sin(lat_rad) | |
lat_cos = math.cos(lat_rad) | |
lat_tan = lat_sin / lat_cos | |
lat_tan2 = lat_tan * lat_tan | |
lat_tan4 = lat_tan2 * lat_tan2 | |
if force_zone_number is None: | |
zone_number = latlon_to_zone_number(latitude, longitude) | |
else: | |
zone_number = force_zone_number | |
zone_letter = latitude_to_zone_letter(latitude) | |
lon_rad = math.radians(longitude) | |
central_lon = zone_number_to_central_longitude(zone_number) | |
central_lon_rad = math.radians(central_lon) | |
n = R / math.sqrt(1 - E * lat_sin**2) | |
c = E_P2 * lat_cos**2 | |
a = lat_cos * (lon_rad - central_lon_rad) | |
a2 = a * a | |
a3 = a2 * a | |
a4 = a3 * a | |
a5 = a4 * a | |
a6 = a5 * a | |
m = R * (M1 * lat_rad - | |
M2 * math.sin(2 * lat_rad) + | |
M3 * math.sin(4 * lat_rad) - | |
M4 * math.sin(6 * lat_rad)) | |
easting = K0 * n * (a + | |
a3 / 6 * (1 - lat_tan2 + c) + | |
a5 / 120 * (5 - 18 * lat_tan2 + lat_tan4 + 72 * c - 58 * E_P2)) + 500000 | |
northing = K0 * (m + n * lat_tan * (a2 / 2 + | |
a4 / 24 * (5 - lat_tan2 + 9 * c + 4 * c**2) + | |
a6 / 720 * (61 - 58 * lat_tan2 + lat_tan4 + 600 * c - 330 * E_P2))) | |
if latitude < 0: | |
northing += 10000000 | |
return easting, northing, zone_number, zone_letter | |
def latitude_to_zone_letter(latitude): | |
for lat_min, zone_letter in ZONE_LETTERS: | |
if latitude >= lat_min: | |
return zone_letter | |
return None | |
def latlon_to_zone_number(latitude, longitude): | |
if 56 <= latitude < 64 and 3 <= longitude < 12: | |
return 32 | |
if 72 <= latitude <= 84 and longitude >= 0: | |
if longitude <= 9: | |
return 31 | |
elif longitude <= 21: | |
return 33 | |
elif longitude <= 33: | |
return 35 | |
elif longitude <= 42: | |
return 37 | |
return int((longitude + 180) / 6) + 1 | |
def zone_number_to_central_longitude(zone_number): | |
return (zone_number - 1) * 6 - 180 + 3 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
almost nothing done here just a simple reminder that tobias and other did a great job