""" utm === .. image:: https://travis-ci.org/Turbo87/utm.png Bidirectional UTM-WGS84 converter for python Usage ----- :: import utm Convert a (latitude, longitude) tuple into an UTM coordinate:: utm.from_latlon(51.2, 7.5) >>> (395201.3103811303, 5673135.241182375, 32, 'U') Convert an UTM coordinate into a (latitude, longitude) tuple:: utm.to_latlon(340000, 5710000, 32, 'U') >>> (51.51852098408468, 6.693872395145327) Speed ----- The library has been compared to the more generic pyproj library by running the unit test suite through pyproj instead of utm. These are the results: * with pyproj (without projection cache): 4.0 - 4.5 sec * with pyproj (with projection cache): 0.9 - 1.0 sec * with utm: 0.4 - 0.5 sec Authors ------- * Tobias Bieniek License ------- Copyright (C) 2012 Tobias Bieniek 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, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ import math __all__ = ['to_latlon', 'from_latlon'] class OutOfRangeError(ValueError): pass 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) _E3 = _E * _E * _E _E4 = _E3 * _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 * _E / 2 - 27 * _E3 / 32) P3 = (21 * _E3 / 16 - 55 * _E4 / 32) P4 = (151 * _E3 / 96) 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): zone_letter = zone_letter.upper() 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 not 'C' <= zone_letter <= 'X' or zone_letter in ['I', 'O']: raise OutOfRangeError('zone letter out of range (must be between C and X)') x = easting - 500000 y = northing if zone_letter < 'N': 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)) 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): 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('northing 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 lon_rad = math.radians(longitude) zone_number = latlon_to_zone_number(latitude, longitude) central_lon = zone_number_to_central_longitude(zone_number) central_lon_rad = math.radians(central_lon) zone_letter = latitude_to_zone_letter(latitude) 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 def haversine(lon1, lat1, lon2, lat2): """ Calculate the great circle distance between two points on the earth (specified in decimal degrees) """ # convert decimal degrees to radians lon1, lat1, lon2, lat2 = map(math.radians, [lon1, lat1, lon2, lat2]) # haversine formula dlon = lon2 - lon1 dlat = lat2 - lat1 a = math.sin(dlat/2)**2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon/2)**2 c = 2 * math.asin(math.sqrt(a)) m = 6367000 * c return m