Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How to convert from UTM to LatLng in python or Javascript

I have a bunch of files with coordinates in UTM form. For each coordinate I have easting, northing and zone. I need to convert this to LatLng for use with Google Map API to show the information in a map.

I have found some online calculators that does this, but no actual code or libraries. http://trac.osgeo.org/proj4js/ is a projection library for Javascript, but looking at the demo it doesn't include UTM projection.

I am still pretty fresh to the entire GIS domain, so what I want is something ala:

(lat,lng) = transform(easting, northing, zone) 
like image 701
Staale Avatar asked Dec 05 '08 13:12

Staale


1 Answers

I ended up finding java code from IBM that solved it: http://www.ibm.com/developerworks/java/library/j-coordconvert/index.html

Just for reference, here is my python implementation of the method I needed:

import math  def utmToLatLng(zone, easting, northing, northernHemisphere=True):     if not northernHemisphere:         northing = 10000000 - northing      a = 6378137     e = 0.081819191     e1sq = 0.006739497     k0 = 0.9996      arc = northing / k0     mu = arc / (a * (1 - math.pow(e, 2) / 4.0 - 3 * math.pow(e, 4) / 64.0 - 5 * math.pow(e, 6) / 256.0))      ei = (1 - math.pow((1 - e * e), (1 / 2.0))) / (1 + math.pow((1 - e * e), (1 / 2.0)))      ca = 3 * ei / 2 - 27 * math.pow(ei, 3) / 32.0      cb = 21 * math.pow(ei, 2) / 16 - 55 * math.pow(ei, 4) / 32     cc = 151 * math.pow(ei, 3) / 96     cd = 1097 * math.pow(ei, 4) / 512     phi1 = mu + ca * math.sin(2 * mu) + cb * math.sin(4 * mu) + cc * math.sin(6 * mu) + cd * math.sin(8 * mu)      n0 = a / math.pow((1 - math.pow((e * math.sin(phi1)), 2)), (1 / 2.0))      r0 = a * (1 - e * e) / math.pow((1 - math.pow((e * math.sin(phi1)), 2)), (3 / 2.0))     fact1 = n0 * math.tan(phi1) / r0      _a1 = 500000 - easting     dd0 = _a1 / (n0 * k0)     fact2 = dd0 * dd0 / 2      t0 = math.pow(math.tan(phi1), 2)     Q0 = e1sq * math.pow(math.cos(phi1), 2)     fact3 = (5 + 3 * t0 + 10 * Q0 - 4 * Q0 * Q0 - 9 * e1sq) * math.pow(dd0, 4) / 24      fact4 = (61 + 90 * t0 + 298 * Q0 + 45 * t0 * t0 - 252 * e1sq - 3 * Q0 * Q0) * math.pow(dd0, 6) / 720      lof1 = _a1 / (n0 * k0)     lof2 = (1 + 2 * t0 + Q0) * math.pow(dd0, 3) / 6.0     lof3 = (5 - 2 * Q0 + 28 * t0 - 3 * math.pow(Q0, 2) + 8 * e1sq + 24 * math.pow(t0, 2)) * math.pow(dd0, 5) / 120     _a2 = (lof1 - lof2 + lof3) / math.cos(phi1)     _a3 = _a2 * 180 / math.pi      latitude = 180 * (phi1 - fact1 * (fact2 + fact3 + fact4)) / math.pi      if not northernHemisphere:         latitude = -latitude      longitude = ((zone > 0) and (6 * zone - 183.0) or 3.0) - _a3      return (latitude, longitude) 

And here I thought it was something simple like easting*x+zone*y or something.

like image 71
Staale Avatar answered Oct 02 '22 02:10

Staale