# Bidirectional UTM-WGS84 converter for Python
# Adapted from the Python script created by Tobias Bieniek
try:
import numpy as mathlib
use_numpy = True
except ImportError:
import math as mathlib
use_numpy = False
__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 = mathlib.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 = "CDEFGHJKLMNPQRSTUVWXX"
def in_bounds(x, lower, upper, upper_strict=False):
if upper_strict and use_numpy:
return lower <= mathlib.min(x) and mathlib.max(x) < upper
elif upper_strict and not use_numpy:
return lower <= x < upper
elif use_numpy:
return lower <= mathlib.min(x) and mathlib.max(x) <= upper
return lower <= x <= upper
def check_valid_zone(zone_number, zone_letter):
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)')
def mixed_signs(x):
return use_numpy and mathlib.min(x) < 0 and mathlib.max(x) >= 0
def negative(x):
if use_numpy:
return mathlib.max(x) < 0
return x < 0
def mod_angle(value):
"""
Returns angle in radians to be between -pi and pi
"""
return (value + mathlib.pi) % (2 * mathlib.pi) - mathlib.pi
[docs]
def to_latlon(easting, northing, zone_number, zone_letter=None, northern=None, strict=True):
"""
This function converts UTM coordinates to Latitude and Longitude
Parameters
----------
easting: int or NumPy array
Easting value of UTM coordinates
northing: int or NumPy array
Northing value of UTM coordinates
zone_number: int
Zone number is represented with global map numbers of a UTM zone
numbers map. For more information `see utmzones <http://www.jaworski.ca/utmzones.htm>`_
zone_letter: str
Zone letter can be represented as string values. UTM zone
designators can be seen in `this reference <http://www.jaworski.ca/utmzones.htm>`_
northern: bool
You can set True or False to set this parameter. Default is None
strict: bool
Raise an OutOfRangeError if outside of bounds
Returns
-------
latitude: float or NumPy array
Latitude between 80 deg S and 84 deg N, e.g. (-80.0 to 84.0)
longitude: float or NumPy array
Longitude between 180 deg W and 180 deg E, e.g. (-180.0 to 180.0).
"""
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 strict:
if not in_bounds(easting, 100000, 1000000, upper_strict=True):
raise OutOfRangeError('easting out of range (must be between 100,000 m and 999,999 m)')
if not in_bounds(northing, 0, 10000000):
raise OutOfRangeError('northing out of range (must be between 0 m and 10,000,000 m)')
check_valid_zone(zone_number, zone_letter)
if zone_letter:
zone_letter = zone_letter.upper()
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 * mathlib.sin(2 * mu) +
P3 * mathlib.sin(4 * mu) +
P4 * mathlib.sin(6 * mu) +
P5 * mathlib.sin(8 * mu))
p_sin = mathlib.sin(p_rad)
p_sin2 = p_sin * p_sin
p_cos = mathlib.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 = mathlib.sqrt(1 - E * p_sin2)
n = R / ep_sin_sqrt
r = (1 - E) / ep_sin
c = E_P2 * 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
longitude = mod_angle(longitude + mathlib.radians(zone_number_to_central_longitude(zone_number)))
return (mathlib.degrees(latitude),
mathlib.degrees(longitude))
[docs]
def from_latlon(latitude, longitude, force_zone_number=None, force_zone_letter=None):
"""
This function converts Latitude and Longitude to UTM coordinate
Parameters
----------
latitude: float or NumPy array
Latitude between 80 deg S and 84 deg N, e.g. (-80.0 to 84.0)
longitude: float or NumPy array
Longitude between 180 deg W and 180 deg E, e.g. (-180.0 to 180.0).
force_zone_number: int
Zone number is represented by global map numbers of an UTM zone
numbers map. You may force conversion to be included within one
UTM zone number. For more information `see utmzones <http://www.jaworski.ca/utmzones.htm>`_
force_zone_letter: str
You may force conversion to be included within one UTM zone
letter. For more information `see utmzones <http://www.jaworski.ca/utmzones.htm>`_
Returns
-------
easting: float or NumPy array
Easting value of UTM coordinates
northing: float or NumPy array
Northing value of UTM coordinates
zone_number: int
Zone number is represented by global map numbers of a UTM zone
numbers map. More information `see utmzones <http://www.jaworski.ca/utmzones.htm>`_
zone_letter: str
Zone letter is represented by a string value. UTM zone designators
can be accessed `here <http://www.jaworski.ca/utmzones.htm>`_
"""
if not in_bounds(latitude, -80.0, 84.0):
raise OutOfRangeError('latitude out of range (must be between 80 deg S and 84 deg N)')
if not in_bounds(longitude, -180.0, 180.0):
raise OutOfRangeError('longitude out of range (must be between 180 deg W and 180 deg E)')
if force_zone_number is not None:
check_valid_zone(force_zone_number, force_zone_letter)
lat_rad = mathlib.radians(latitude)
lat_sin = mathlib.sin(lat_rad)
lat_cos = mathlib.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
if force_zone_letter is None:
zone_letter = latitude_to_zone_letter(latitude)
else:
zone_letter = force_zone_letter
lon_rad = mathlib.radians(longitude)
central_lon = zone_number_to_central_longitude(zone_number)
central_lon_rad = mathlib.radians(central_lon)
n = R / mathlib.sqrt(1 - E * lat_sin**2)
c = E_P2 * lat_cos**2
a = lat_cos * mod_angle(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 * mathlib.sin(2 * lat_rad) +
M3 * mathlib.sin(4 * lat_rad) -
M4 * mathlib.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 mixed_signs(latitude):
raise ValueError("latitudes must all have the same sign")
elif negative(latitude):
northing += 10000000
return easting, northing, zone_number, zone_letter
def latitude_to_zone_letter(latitude):
# If the input is a numpy array, just use the first element
# User responsibility to make sure that all points are in one zone
if use_numpy and isinstance(latitude, mathlib.ndarray):
latitude = latitude.flat[0]
if -80 <= latitude <= 84:
return ZONE_LETTERS[int(latitude + 80) >> 3]
else:
return None
def latlon_to_zone_number(latitude, longitude):
# If the input is a numpy array, just use the first element
# User responsibility to make sure that all points are in one zone
if use_numpy:
if isinstance(latitude, mathlib.ndarray):
latitude = latitude.flat[0]
if isinstance(longitude, mathlib.ndarray):
longitude = longitude.flat[0]
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
class OutOfRangeError(ValueError):
pass