How to convert from UTM to LatLng in python or Jav

问题:

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)

回答1:

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.

回答2:

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.

回答3:

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.

回答4:

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/[e-mail]/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"])

回答5:

You could use Proj4js, as follows.

Download Proj4JS from GitHub, using this link.

The following code will convert from UTM to longitude latitude

<html>
<head>
  <script src="proj4.js"></script>

  <script>
    var utm = "+proj=utm +zone=32";
    var wgs84 = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs";
    console.log(proj4(utm,wgs84,[539884, 4942158]));
  </script>
</head>
<body>

</body>
</html>

In this code, the UTM zone is 32, as should be obvious. The Easting is 539884, and the Northing is 4942158. The result is:

[9.502832656648073, 44.631671014204365] 

Which is to say 44.631671014204365N, 9.502832656648073E. Which I have verified is correct.

If you need other projections, you can find their strings here.

回答6:

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;
      }

回答7:

////////////////////////////////////////////////////////////////////////////////////////////
//
// ToLL - function to compute Latitude and Longitude given UTM Northing and Easting in meters
//
//  Description:
//    This member function converts input north and east coordinates
//    to the corresponding Northing and Easting values relative to the defined
//    UTM zone.  Refer to the reference in this file's header.
//
//  Parameters:
//    north   - (i) Northing (meters)
//    east    - (i) Easting (meters)
//    utmZone - (i) UTM Zone of the North and East parameters
//    lat     - (o) Latitude in degrees 
//    lon     - (o) Longitude in degrees
//
function ToLL(north,east,utmZone)
{ 
  // This is the lambda knot value in the reference
  var LngOrigin = DegToRad(utmZone * 6 - 183)

  // The following set of class constants define characteristics of the
  // ellipsoid, as defined my the WGS84 datum.  These values need to be
  // changed if a different dataum is used.    

  var FalseNorth = 0.  // South or North?
  //if (lat < 0.) FalseNorth = 10000000.  // South or North?
  //else          FalseNorth = 0.   

  var Ecc = 0.081819190842622       // Eccentricity
  var EccSq = Ecc * Ecc
  var Ecc2Sq = EccSq / (1. - EccSq)
  var Ecc2 = Math.sqrt(Ecc2Sq)      // Secondary eccentricity
  var E1 = ( 1 - Math.sqrt(1-EccSq) ) / ( 1 + Math.sqrt(1-EccSq) )
  var E12 = E1 * E1
  var E13 = E12 * E1
  var E14 = E13 * E1

  var SemiMajor = 6378137.0         // Ellipsoidal semi-major axis (Meters)
  var FalseEast = 500000.0          // UTM East bias (Meters)
  var ScaleFactor = 0.9996          // Scale at natural origin

  // Calculate the Cassini projection parameters

  var M1 = (north - FalseNorth) / ScaleFactor
  var Mu1 = M1 / ( SemiMajor * (1 - EccSq/4.0 - 3.0*EccSq*EccSq/64.0 -
    5.0*EccSq*EccSq*EccSq/256.0) )

  var Phi1 = Mu1 + (3.0*E1/2.0 - 27.0*E13/32.0) * Math.sin(2.0*Mu1)
    + (21.0*E12/16.0 - 55.0*E14/32.0)           * Math.sin(4.0*Mu1)
    + (151.0*E13/96.0)                          * Math.sin(6.0*Mu1)
    + (1097.0*E14/512.0)                        * Math.sin(8.0*Mu1)

  var sin2phi1 = Math.sin(Phi1) * Math.sin(Phi1)
  var Rho1 = (SemiMajor * (1.0-EccSq) ) / Math.pow(1.0-EccSq*sin2phi1,1.5)
  var Nu1 = SemiMajor / Math.sqrt(1.0-EccSq*sin2phi1)

  // Compute parameters as defined in the POSC specification.  T, C and D

  var T1 = Math.tan(Phi1) * Math.tan(Phi1)
  var T12 = T1 * T1
  var C1 = Ecc2Sq * Math.cos(Phi1) * Math.cos(Phi1)
  var C12 = C1 * C1
  var D  = (east - FalseEast) / (ScaleFactor * Nu1)
  var D2 = D * D
  var D3 = D2 * D
  var D4 = D3 * D
  var D5 = D4 * D
  var D6 = D5 * D

  // Compute the Latitude and Longitude and convert to degrees
  var lat = Phi1 - Nu1*Math.tan(Phi1)/Rho1 *
    ( D2/2.0 - (5.0 + 3.0*T1 + 10.0*C1 - 4.0*C12 - 9.0*Ecc2Sq)*D4/24.0
     + (61.0 + 90.0*T1 + 298.0*C1 + 45.0*T12 - 252.0*Ecc2Sq - 3.0*C12)*D6/720.0 )

  lat = RadToDeg(lat)

  var lon = LngOrigin + 
    ( D - (1.0 + 2.0*T1 + C1)*D3/6.0
      + (5.0 - 2.0*C1 + 28.0*T1 - 3.0*C12 + 8.0*Ecc2Sq + 24.0*T12)*D5/120.0) / Math.cos(Phi1)

  lon = RadToDeg(lon)

  // Create a object to store the calculated Latitude and Longitude values
  var sendLatLon = new PC_LatLon(lat,lon)

  // Returns a PC_LatLon object
  return sendLatLon
}

////////////////////////////////////////////////////////////////////////////////////////////
//
//  RadToDeg - function that inputs a value in radians and returns a value in degrees
//
function RadToDeg(value)
{
  return ( value * 180.0 / Math.PI )
}

////////////////////////////////////////////////////////////////////////////////////////////
//
// PC_LatLon - this psuedo class is used to store lat/lon values computed by the ToLL 
//  function.
//
function PC_LatLon(inLat,inLon)
{
  this.lat       = inLat     // Store Latitude in decimal degrees
  this.lon       = inLon     // Store Longitude in decimal degrees
}

回答8:

There is a perl module via CPAN called Geography::NationalGrid which can convert easting/northing to lat/longs. That may help.

Alternatively there are lots of scripts on the movable-type site that let you convert lat/long and easting/northings.

回答9:

One problem I had with using proj4js was that it needed the exact zone as @Richard points out. I found a great resource here which can convert WGS to UTM and wrote a cleaner wrapper in JavaScript:

https://github.com/urbanetic/utm-converter

回答10:

The answer by Staale worked for me with a small modification - The math module cannot handle Pandas Series, so I replaced all math functions with numpy.

However, checking in QGIS, I see about 4m difference between the UTM and LAT/LON coordinates.

Code below:

import numpy as np

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 - np.power(e, 2) / 4.0 - 3 * np.power(e, 4) / 64.0 - 5 * np.power(e, 6) / 256.0))

ei = (1 - np.power((1 - e * e), (1 / 2.0))) / (1 + np.power((1 - e * e), (1 / 2.0)))

ca = 3 * ei / 2 - 27 * np.power(ei, 3) / 32.0

cb = 21 * np.power(ei, 2) / 16 - 55 * np.power(ei, 4) / 32
cc = 151 * np.power(ei, 3) / 96
cd = 1097 * np.power(ei, 4) / 512
phi1 = mu + ca * np.sin(2 * mu) + cb * np.sin(4 * mu) + cc * np.sin(6 * mu) + cd * np.sin(8 * mu)

n0 = a / np.power((1 - np.power((e * np.sin(phi1)), 2)), (1 / 2.0))

r0 = a * (1 - e * e) / np.power((1 - np.power((e * np.sin(phi1)), 2)), (3 / 2.0))
fact1 = n0 * np.tan(phi1) / r0

_a1 = 500000 - easting
dd0 = _a1 / (n0 * k0)
fact2 = dd0 * dd0 / 2

t0 = np.power(np.tan(phi1), 2)
Q0 = e1sq * np.power(np.cos(phi1), 2)
fact3 = (5 + 3 * t0 + 10 * Q0 - 4 * Q0 * Q0 - 9 * e1sq) * np.power(dd0, 4) / 24

fact4 = (61 + 90 * t0 + 298 * Q0 + 45 * t0 * t0 - 252 * e1sq - 3 * Q0 * Q0) * np.power(dd0, 6) / 720

lof1 = _a1 / (n0 * k0)
lof2 = (1 + 2 * t0 + Q0) * np.power(dd0, 3) / 6.0
lof3 = (5 - 2 * Q0 + 28 * t0 - 3 * np.power(Q0, 2) + 8 * e1sq + 24 * np.power(t0, 2)) * np.power(dd0, 5) / 120
_a2 = (lof1 - lof2 + lof3) / np.cos(phi1)
_a3 = _a2 * 180 / np.pi

latitude = 180 * (phi1 - fact1 * (fact2 + fact3 + fact4)) / np.pi

if not northernHemisphere:
    latitude = -latitude

longitude = ((zone > 0) and (6 * zone - 183.0) or 3.0) - _a3

return (latitude, longitude)

That way I can do this directly:

df['LAT'], df['LON']=utmToLatLng(31, df['X'], df['Y'], northernHemisphere=True)