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)
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.
What I found is the following site: http://home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html It has a javascript converter, you should check the algorithm there. From the page:
Programmers: The JavaScript source code in this document may be copied and reused without restriction.
According to this page, UTM is supported by proj4js.
http://trac.osgeo.org/proj4js/wiki/UserGuide#Supportedprojectionclasses
You may also want to take a look at GDAL. The gdal library has excellent python support, though it may be a bit overkill if you're only doing projection conversion.
I'm new to this as well and have been studying up on the subject recently.
Here's a method I found using the python gdal pacakge (the osr package is included in gdal). The gdal package is pretty powerful, but the documentation could be better.
This is derived from a discussion here: http://www.mail-archive.com/[email protected]/msg12398.html
import osr
def transform_utm_to_wgs84(easting, northing, zone):
utm_coordinate_system = osr.SpatialReference()
utm_coordinate_system.SetWellKnownGeogCS("WGS84") # Set geographic coordinate system to handle lat/lon
is_northern = northing > 0
utm_coordinate_system.SetUTM(zone, is_northern)
wgs84_coordinate_system = utm_coordinate_system.CloneGeogCS() # Clone ONLY the geographic coordinate system
# create transform component
utm_to_wgs84_transform = osr.CoordinateTransformation(utm_coordinate_system, wgs84_coordinate_system) # (<from>, <to>)
return utm_to_wgs84_transform.TransformPoint(easting, northing, 0) # returns lon, lat, altitude
And here's the method for converting from a lat, lon in wgs84 (what most gps units report) to utm:
def transform_wgs84_to_utm(lon, lat):
def get_utm_zone(longitude):
return (int(1+(longitude+180.0)/6.0))
def is_northern(latitude):
"""
Determines if given latitude is a northern for UTM
"""
if (latitude < 0.0):
return 0
else:
return 1
utm_coordinate_system = osr.SpatialReference()
utm_coordinate_system.SetWellKnownGeogCS("WGS84") # Set geographic coordinate system to handle lat/lon
utm_coordinate_system.SetUTM(get_utm_zone(lon), is_northern(lat))
wgs84_coordinate_system = utm_coordinate_system.CloneGeogCS() # Clone ONLY the geographic coordinate system
# create transform component
wgs84_to_utm_transform = osr.CoordinateTransformation(wgs84_coordinate_system, utm_coordinate_system) # (<from>, <to>)
return wgs84_to_utm_transform.TransformPoint(lon, lat, 0) # returns easting, northing, altitude
I also found that if you've already got django/gdal installed and you know the EPSG code for the UTM zone you're working on, you can just use the Point()
transform() method.
from django.contrib.gis.geos import Point
utm2epsg = {"54N": 3185, ...}
p = Point(lon, lat, srid=4326) # 4326 = WGS84 epsg code
p.transform(utm2epsg["54N"])
A Javascript version of Staale answer
function utmToLatLng(zone, easting, northing, northernHemisphere){
if (!northernHemisphere){
northing = 10000000 - northing;
}
var a = 6378137;
var e = 0.081819191;
var e1sq = 0.006739497;
var k0 = 0.9996;
var arc = northing / k0;
var mu = arc / (a * (1 - Math.pow(e, 2) / 4.0 - 3 * Math.pow(e, 4) / 64.0 - 5 * Math.pow(e, 6) / 256.0));
var ei = (1 - Math.pow((1 - e * e), (1 / 2.0))) / (1 + Math.pow((1 - e * e), (1 / 2.0)));
var ca = 3 * ei / 2 - 27 * Math.pow(ei, 3) / 32.0;
var cb = 21 * Math.pow(ei, 2) / 16 - 55 * Math.pow(ei, 4) / 32;
var cc = 151 * Math.pow(ei, 3) / 96;
var cd = 1097 * Math.pow(ei, 4) / 512;
var phi1 = mu + ca * Math.sin(2 * mu) + cb * Math.sin(4 * mu) + cc * Math.sin(6 * mu) + cd * Math.sin(8 * mu);
var n0 = a / Math.pow((1 - Math.pow((e * Math.sin(phi1)), 2)), (1 / 2.0));
var r0 = a * (1 - e * e) / Math.pow((1 - Math.pow((e * Math.sin(phi1)), 2)), (3 / 2.0));
var fact1 = n0 * Math.tan(phi1) / r0;
var _a1 = 500000 - easting;
var dd0 = _a1 / (n0 * k0);
var fact2 = dd0 * dd0 / 2;
var t0 = Math.pow(Math.tan(phi1), 2);
var Q0 = e1sq * Math.pow(Math.cos(phi1), 2);
var fact3 = (5 + 3 * t0 + 10 * Q0 - 4 * Q0 * Q0 - 9 * e1sq) * Math.pow(dd0, 4) / 24;
var fact4 = (61 + 90 * t0 + 298 * Q0 + 45 * t0 * t0 - 252 * e1sq - 3 * Q0 * Q0) * Math.pow(dd0, 6) / 720;
var lof1 = _a1 / (n0 * k0);
var lof2 = (1 + 2 * t0 + Q0) * Math.pow(dd0, 3) / 6.0;
var lof3 = (5 - 2 * Q0 + 28 * t0 - 3 * Math.pow(Q0, 2) + 8 * e1sq + 24 * Math.pow(t0, 2)) * Math.pow(dd0, 5) / 120;
var _a2 = (lof1 - lof2 + lof3) / Math.cos(phi1);
var _a3 = _a2 * 180 / Math.PI;
var latitude = 180 * (phi1 - fact1 * (fact2 + fact3 + fact4)) / Math.PI;
if (!northernHemisphere){
latitude = -latitude;
}
var longitude = ((zone > 0) && (6 * zone - 183.0) || 3.0) - _a3;
var obj = {
latitude : latitude,
longitude: longitude
};
return obj;
}