mosp.geo.utm

1 # -*- coding: utf-8 -*- 2 """Working with UTM and Lat/Long coordinates 3 4 Based on http://home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html 5 6 """ 7 8 from math import sin, cos, sqrt, tan, pi, floor 9 10 __author__ = "P. Tute, C. Taylor" 11 __maintainer__ = "B. Henne" 12 __contact__ = "henne@dcsec.uni-hannover.de" 13 __copyright__ = "(c) 1997-2003 C. Taylor, 2010-2011, DCSec, Leibniz Universitaet Hannover, Germany" 14 __license__ = """ 15 The Python code is based on Javascript code of Charles L. Taylor from 16 http://home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html 17 The author allowed us to reuse his code. 18 Original license: 19 The documents, images, and other materials on the Chuck Taylor 20 Web site are copyright (c) 1997-2003 by C. Taylor, unless otherwise noted. 21 Visitors are permitted to download the materials located on the Chuck Taylor 22 Web site to their own systems, on a temporary basis, for the purpose of viewing 23 them or, in the case of software applications, using them in accordance with 24 their intended purpose. 25 Republishing or redistribution of any of the materials located on the Chuck Taylor 26 Web site requires permission from C. Taylor. MoSP has been granted to use the code.""" 27 28 29 # Ellipsoid model constants (actual values here are for WGS84) 30 UTMScaleFactor = 0.9996 #: Ellipsoid model constant 31 sm_a = 6378137.0 #: Ellipsoid model constant: Semi-major axis a 32 sm_b = 6356752.314 #: Ellipsoid model constant: Semi-major axis b 33 sm_EccSquared = 6.69437999013e-03 #: Ellipsoid model constant 34 35

36 -def long_to_zone(lon):

37 """Calculates the current UTM-zone for a given longitude.""" 38 return floor((lon + 180.0) / 6) + 1

39 40

41 -def rad_to_deg(rad):

42 """Converts radians to degrees.""" 43 return (rad / pi * 180.0)

44 45

46 -def deg_to_rad(deg):

47 """Converts degrees to radians.""" 48 return (deg / 180.0 * pi)

49 50

51 -def UTMCentralMeridian(zone):

52 """Calculates the central meridian for the given UTM-zone.""" 53 return deg_to_rad(-183.0 + (zone * 6.0))

54 55

56 -def ArcLengthOfMeridian(phi):

57 """Computes the ellipsoidal distance from the equator to a point at a 58 given latitude. 59 60 Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J., 61 GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994.""" 62 63 # Precalculate n 64 n = (sm_a - sm_b) / (sm_a + sm_b) 65 66 # Precalculate alpha 67 alpha = (((sm_a + sm_b) / 2.0) 68 * (1.0 + (n**2.0 / 4.0) + (n**4.0 / 64.0))) 69 70 # Precalculate beta 71 beta = ((-3.0 * n / 2.0) + (9.0 * n**3.0 / 16.0) 72 + (-3.0 * n**5.0 / 32.0)) 73 74 # Precalculate gamma 75 gamma = ((15.0 * n**2.0 / 16.0) 76 + (-15.0 * n**4.0 / 32.0)) 77 78 # Precalculate delta 79 delta = ((-35.0 * n**3.0 / 48.0) 80 + (105.0 * n**5.0 / 256.0)) 81 82 # Precalculate epsilon 83 epsilon = (315.0 * n**4.0 / 512.0) 84 85 # Now calculate the sum of the series and return 86 result = (alpha 87 * (phi + (beta * sin(2.0 * phi)) 88 + (gamma * sin(4.0 * phi)) 89 + (delta * sin(6.0 * phi)) 90 + (epsilon * sin(8.0 * phi)))) 91 92 return result

93 94

95 -def MapLatLonToXY(phi, lambd, lambd0, xy):

96 """Converts a latitude/longitude pair to x and y coordinates in the 97 Transverse Mercator projection. Note that Transverse Mercator is not 98 the same as UTM a scale factor is required to convert between them. 99 100 Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J., 101 GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994.""" 102 103 # Precalculate ep2 104 ep2 = (sm_a**2.0 - sm_b**2.0) / sm_b**2.0 105 106 # Precalculate nu2 107 nu2 = ep2 * cos(phi)**2.0 108 109 # Precalculate N 110 N = sm_a**2.0 / (sm_b * sqrt(1 + nu2)) 111 112 # Precalculate t 113 t = tan(phi) 114 t2 = t**2.0 115 116 # Precalculate l 117 l = lambd - lambd0 118 119 # Precalculate coefficients for l**n in the equations below 120 # so a normal human being can read the expressions for easting 121 # and northing 122 # -- l**1 and l**2 have coefficients of 1.0 123 l3coef = 1.0 - t2 + nu2 124 125 l4coef = 5.0 - t2 + 9 * nu2 + 4.0 * (nu2**2.0) 126 127 l5coef = 5.0 - 18.0 * t2 + (t2**2.0) + 14.0 * nu2 - 58.0 * t2 * nu2 128 129 l6coef = 61.0 - 58.0 * t2 + (t2**2.0) + 270.0 * nu2 - 330.0 * t2 * nu2 130 131 l7coef = 61.0 - 479.0 * t2 + 179.0 * (t2**2.0) - (t2**3.0) 132 133 l8coef = 1385.0 - 3111.0 * t2 + 543.0 * (t2**2.0) - (t2**3.0) 134 135 # Calculate easting (x) 136 xy[0] = (N * cos(phi) * l + (N / 6.0 * cos(phi)**3.0 * l3coef * l**3.0) 137 + (N / 120.0 * cos(phi)**5.0 * l5coef * l**5.0) 138 + (N / 5040.0 * cos(phi)**7.0 * l7coef * l**7.0)) 139 140 # Calculate northing (y) 141 xy[1] = (ArcLengthOfMeridian(phi) 142 + (t / 2.0 * N * cos(phi)**2.0 * l**2.0) 143 + (t / 24.0 * N * cos(phi)**4.0 * l4coef * l**4.0) 144 + (t / 720.0 * N * cos(phi)**6.0 * l6coef * l**6.0) 145 + (t / 40320.0 * N * cos(phi)**8.0 * l8coef * l**8.0)) 146 147 return

148 149

150 -def latlong_to_utm(lon, lat, zone = None):

151 """Converts a latitude/longitude pair to x and y coordinates in the 152 Universal Transverse Mercator projection.""" 153 154 if zone is None: 155 zone = long_to_zone(lon) 156 157 xy = [0, 0] 158 MapLatLonToXY(deg_to_rad(lat), deg_to_rad(lon), UTMCentralMeridian(zone), xy) 159 160 xy[0] = xy[0] * UTMScaleFactor + 500000.0 161 xy[1] = xy[1] * UTMScaleFactor 162 163 if xy[1] < 0.0: 164 xy[1] = xy[1] + 10000000.0 165 166 return [round(coord, 2) for coord in xy]

167 168

169 -def FootpointLatitude(y):

170 """Computes the footpoint latitude for use in converting transverse 171 Mercator coordinates to ellipsoidal coordinates. 172 173 Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J., 174 GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994.""" 175 176 # Precalculate n (Eq. 10.18) 177 n = (sm_a - sm_b) / (sm_a + sm_b) 178 179 # Precalculate alpha_ (Eq. 10.22) 180 # (Same as alpha in Eq. 10.17) 181 alpha_ = (((sm_a + sm_b) / 2.0) 182 * (1 + (n**2.0 / 4) + (n**4.0 / 64))) 183 184 # Precalculate y_ (Eq. 10.23) 185 y_ = y / alpha_ 186 187 # Precalculate beta_ (Eq. 10.22) 188 beta_ = ((3.0 * n / 2.0) + (-27.0 * n**3.0 / 32.0) 189 + (269.0 * n**5.0 / 512.0)) 190 191 # Precalculate gamma_ (Eq. 10.22) 192 gamma_ = ((21.0 * n**2.0 / 16.0) 193 + (-55.0 * n**4.0 / 32.0)) 194 195 # Precalculate delta_ (Eq. 10.22) 196 delta_ = ((151.0 * n**3.0 / 96.0) 197 + (-417.0 * n**5.0 / 128.0)) 198 199 # Precalculate epsilon_ (Eq. 10.22) 200 epsilon_ = ((1097.0 * n**4.0 / 512.0)) 201 202 # Now calculate the sum of the series (Eq. 10.21) 203 result = (y_ + (beta_ * sin(2.0 * y_)) 204 + (gamma_ * sin(4.0 * y_)) 205 + (delta_ * sin(6.0 * y_)) 206 + (epsilon_ * sin(8.0 * y_))) 207 208 return result

209 210

211 -def MapXYToLatLon(x, y, lambda0):

212 """Converts x and y coordinates in the Transverse Mercator projection to 213 a latitude/longitude pair. Note that Transverse Mercator is not 214 the same as UTM a scale factor is required to convert between them. 215 216 Reference: Hoffmann-Wellenhof, B., Lichtenegger, H., and Collins, J., 217 GPS: Theory and Practice, 3rd ed. New York: Springer-Verlag Wien, 1994. 218 219 Remarks: 220 The local variables Nf, nuf2, tf, and tf2 serve the same purpose as 221 N, nu2, t, and t2 in MapLatLonToXY, but they are computed with respect 222 sto the footpoint latitude phif. 223 224 x1frac, x2frac, x2poly, x3poly, etc. are to enhance readability and 225 to optimize computations.""" 226 227 philambda = [] 228 229 # Get the value of phif, the footpoint latitude. 230 phif = FootpointLatitude(y) 231 232 # Precalculate ep2 233 ep2 = ((sm_a**2.0 - sm_b**2.0) 234 / sm_b**2.0) 235 236 # Precalculate cos (phif) 237 cf = cos(phif) 238 239 # Precalculate nuf2 240 nuf2 = ep2 * cf**2.0 241 242 # Precalculate Nf and initialize Nfpow 243 Nf = sm_a**2.0 / (sm_b * sqrt(1 + nuf2)) 244 Nfpow = Nf 245 246 # Precalculate tf 247 tf = tan(phif) 248 tf2 = tf**2 249 tf4 = tf2**2 250 251 # Precalculate fractional coefficients for x**n in the equations 252 # below to simplify the expressions for latitude and longitude. 253 x1frac = 1.0 / (Nfpow * cf) 254 255 Nfpow *= Nf # now equals Nf**2) 256 x2frac = tf / (2.0 * Nfpow) 257 258 Nfpow *= Nf # now equals Nf**3) 259 x3frac = 1.0 / (6.0 * Nfpow * cf) 260 261 Nfpow *= Nf # now equals Nf**4) 262 x4frac = tf / (24.0 * Nfpow) 263 264 Nfpow *= Nf # now equals Nf**5) 265 x5frac = 1.0 / (120.0 * Nfpow * cf) 266 267 Nfpow *= Nf # now equals Nf**6) 268 x6frac = tf / (720.0 * Nfpow) 269 270 Nfpow *= Nf # now equals Nf**7) 271 x7frac = 1.0 / (5040.0 * Nfpow * cf) 272 273 Nfpow *= Nf # now equals Nf**8) 274 x8frac = tf / (40320.0 * Nfpow) 275 276 # Precalculate polynomial coefficients for x**n. 277 # -- x**1 does not have a polynomial coefficient. 278 x2poly = -1.0 - nuf2 279 280 x3poly = -1.0 - 2 * tf2 - nuf2 281 282 x4poly = (5.0 + 3.0 * tf2 + 6.0 * nuf2 - 6.0 * tf2 * nuf2 283 - 3.0 * (nuf2 *nuf2) - 9.0 * tf2 * (nuf2 * nuf2)) 284 285 x5poly = 5.0 + 28.0 * tf2 + 24.0 * tf4 + 6.0 * nuf2 + 8.0 * tf2 * nuf2 286 287 x6poly = (-61.0 - 90.0 * tf2 - 45.0 * tf4 - 107.0 * nuf2 288 + 162.0 * tf2 * nuf2) 289 290 x7poly = -61.0 - 662.0 * tf2 - 1320.0 * tf4 - 720.0 * (tf4 * tf2) 291 292 x8poly = 1385.0 + 3633.0 * tf2 + 4095.0 * tf4 + 1575 * (tf4 * tf2) 293 294 # Calculate latitude 295 philambda.append(phif + x2frac * x2poly * x**2 296 + x4frac * x4poly * x**4.0 297 + x6frac * x6poly * x**6.0 298 + x8frac * x8poly * x**8.0) 299 300 # Calculate longitude 301 philambda.append(lambda0 + x1frac * x 302 + x3frac * x3poly * x**3.0 303 + x5frac * x5poly * x**5.0 304 + x7frac * x7poly * x**7.0) 305 306 return philambda

307 308

309 -def utm_to_latlong(x, y, zone, southhemi=False):

310 """Converts x and y coordinates in the Universal Transverse Mercator 311 projection to a latitude/longitude pair.""" 312 313 x -= 500000.0 314 x /= UTMScaleFactor 315 316 # If in southern hemisphere, adjust y accordingly. 317 if southhemi: 318 y -= 10000000.0 319 320 y /= UTMScaleFactor 321 322 cmeridian = UTMCentralMeridian(zone) 323 324 lat_lon = MapXYToLatLon(x, y, cmeridian) 325 return list(reversed([rad_to_deg(i) for i in lat_lon]))

326

Source Code for Module mosp.geo.utm