如何在 python 或 Javascript 中将 UTM 转换为 LatLng
Posted
技术标签:
【中文标题】如何在 python 或 Javascript 中将 UTM 转换为 LatLng【英文标题】：How to convert from UTM to LatLng in python or Javascript 【发布时间】：2010-09-25 12:35:19 【问题描述】：我有一堆带有 UTM 形式坐标的文件。对于每个坐标，我都有东、北和区域。我需要将其转换为 LatLng 以与 Google Map API 一起使用以在地图中显示信息。
我找到了一些可以执行此操作的在线计算器，但没有实际的代码或库。 http://trac.osgeo.org/proj4js/ 是 javascript 的投影库，但看演示它不包括 UTM 投影。
我对整个 GIS 领域还是很陌生，所以我想要的是 ala：
(lat,lng) = transform(easting, northing, zone)
【问题讨论】：
【参考方案1】：我最终从 IBM 找到了解决该问题的 java 代码：http://www.ibm.com/developerworks/java/library/j-coordconvert/index.html
仅供参考，这里是我需要的方法的python实现：
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)
在这里，我认为它很简单，例如 easting*x+zone*y 之类的。
【讨论】：
我不怪你认为这很简单。我上周从纬度/经度转换为米勒投影的东/北，我仍在里面庆祝。 为什么你的函数需要区域？你不应该只需要东、北、半球吗？ @Franki：因为 UTM 坐标是相对于预定义区域的。 UTM on Wikipedia. 你拯救了我的一天，我一直在寻找这样的东西将近 2 天，而且效果很好，我做了一个你的脚本的 javascript 版本，如果其他人觉得它有帮助，我会作为答案...非常感谢您【参考方案2】：我找到的是以下网站：http://home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html 它有一个javascript转换器，你应该在那里检查算法。从页面：
程序员：本文档中的 JavaScript 源代码可以不受限制地复制和重复使用。
【讨论】：
这是我见过的最干净的 Javascript 版本。 正如@monkut 的回答中提到的，python.org 包括来自 python 的 GDAL 的 Swig 绑定：pypi.python.org/pypi/GDAL 链接好像坏了！【参考方案3】：根据this page，proj4js支持UTM。
http://trac.osgeo.org/proj4js/wiki/UserGuide#Supportedprojectionclasses
您可能还想看看GDAL。 gdal 库具有出色的 python 支持，但如果您只进行投影转换，它可能有点矫枉过正。
【讨论】：
+1 PROJ4 几乎可以做任何你梦寐以求的事情，所以如果 proj4js 是一个真正的移植版，它也可以做到。 +1 proj4s 是要走的路。没有必要重新发明***.. proj4s 从可以添加任何投影的配置文件中读取 - 将引用键入spatialreference.org 以获取 proj4js 字符串，例如spatialreference.org/ref/epsg/4326/proj4js【参考方案4】：我也是这方面的新手，最近一直在研究这个主题。
这是我使用 python gdal pacakge 找到的一种方法（osr 包包含在 gdal 中）。 gdal 包非常强大，但文档可能会更好。
这是从这里的讨论中得出的： http://www.mail-archive.com/gdal-dev@lists.osgeo.org/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
以下是从 wgs84（大多数 gps 单位报告的）中的 lat、lon 转换为 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
我还发现，如果您已经安装了 django/gdal，并且知道您正在处理的 UTM 区域的 EPSG 代码，则可以使用 Point() transform() 方法。
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】：Staale 答案的 Javascript 版本
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;
【讨论】：
【参考方案6】：您可以使用 Proj4js，如下所示。
从 GitHub 下载 Proj4JS，使用 this 链接。
以下代码将从UTM转换为经度纬度
<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>
在这段代码中，UTM 区域是 32，这应该很明显。东距为 539884，北距为 4942158。结果为：
[9.502832656648073, 44.631671014204365]
也就是说 44.631671014204365N, 9.502832656648073E。我的verified 是正确的。
如果你需要其他的投影，你可以找到他们的字符串here。
【讨论】：
赤道以南与否【参考方案7】：有一个Python库（utm）“Bidirectional UTM-WGS84 converter for python”，可以管理单个坐标点或系列。
将（纬度，经度）元组转换为UTM坐标：
>>> utm.from_latlon(51.2, 7.5)
结果
(395201.3103811303, 5673135.241182375, 32, 'U')
返回的格式为(EASTING, NORTHING, ZONE_NUMBER, ZONE_LETTER)。
如果使用系列，结果 EASTING 和 NORTHING 将具有相同的形状。
>>> utm.from_latlon(np.array([51.2, 49.0]), np.array([7.5, 8.4]))
结果
(array([395201.31038113, 456114.59586214]), array([5673135.24118237, 5427629.20426126]), 32, 'U')
反之，将 UTM 坐标转换为（纬度，经度）元组：
>>> utm.to_latlon(340000, 5710000, 32, 'U')
(51.51852098408468, 6.693872395145327)
使用语法utm.to_latlon(EASTING, NORTHING, ZONE_NUMBER, ZONE_LETTER)。如果要使用系列，在这种情况下也是如此。
根据文档，传递一系列比多次单点调用函数更快。
【讨论】：
【参考方案8】：////////////////////////////////////////////////////////////////////////////////////////////
//
// 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
【讨论】：
DegToRad 方法丢失【参考方案9】：我在使用 proj4js 时遇到的一个问题是它需要 @Richard 指出的确切区域。我找到了一个很棒的资源here，它可以将 WGS 转换为 UTM，并用 JavaScript 编写了一个更干净的包装器：
https://github.com/urbanetic/utm-converter
【讨论】：
【参考方案10】：通过 CPAN 有一个名为 Geography::NationalGrid 的 perl 模块，它可以将东/北转换为纬度/经度。这可能会有所帮助。
另外，movable-type site 上有很多脚本可以让您转换纬度/经度和东向/北向。
【讨论】：
【参考方案11】：Staale 的回答对我进行了小修改 - 数学模块无法处理 Pandas 系列，所以我用 numpy 替换了所有数学函数。
但是，检查 QGIS，我发现 UTM 和 LAT/LON 坐标之间大约有 4m 的差异。
代码如下：
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)
这样我可以直接这样做：
df['LAT'], df['LON']=utmToLatLng(31, df['X'], df['Y'], northernHemisphere=True)
【讨论】：
【参考方案12】：将@TreyA 转换为 kotlin（Java 库中的一些数学函数），仅针对它的情况。似乎比@Staale 和@sandino 更准确一些
fun utmToLatLng2(
north:Double?,
east:Double?,
zone:Long,
northernHemisphere: Boolean
):LatLng?
if((north==null)||(east == null)) return null
val lngOrigin = Math.toRadians(zone*6.0 - 183.0)
val falseNorth = if(northernHemisphere) 0.toDouble() else 10000000.toDouble()
val ecc = 0.081819190842622
val eccSq = ecc*ecc
val ecc2Sq = eccSq / (1.0 - eccSq)
//var ecc2 = sqrt(ecc2Sq) //not in use ?
val e1 = (1.0 - sqrt(1.0-eccSq))/(1.0 + sqrt(1-eccSq))
val e12 = e1*e1
val e13 = e12*e1
val e14 = e13*e1
val semiMajor = 6378137.0
val falseEast = 500000.0
val scaleFactor = 0.9996
//Cassini
val m1 = (north - falseNorth) / scaleFactor
val mu1 = m1 / (semiMajor * (1.0 - eccSq/4.0 - 3.0*eccSq*eccSq/64.0 - 5.0*eccSq*eccSq*eccSq/256.0))
val phi1 = mu1 +
(3.0 * e1 / 2.0 - 27.0 * e13 / 32.0) * sin(2.0 * mu1) +
(21.0 * e12 / 16.0 - 55.0 * e14 / 32.0) * sin(4.0 * mu1) +
(151.0 * e13 / 96.0) * sin(6.0 * mu1) +
(1097.0 * e14 / 512.0) * sin( 8.0 * mu1)
val sin2phi1 = sin(phi1) * sin(phi1)
val rho1 = (semiMajor * (1.0-eccSq)) / (1.0 - eccSq * sin2phi1).pow(1.5)
val nu1 = semiMajor / sqrt(1.0-eccSq*sin2phi1)
//POSC
val t1 = tan(phi1)*tan(phi1)
val t12 = t1*t1
val c1 = ecc2Sq*cos(phi1)*cos(phi1)
val c12 = c1*c1
val d = (east - falseEast) / (scaleFactor*nu1)
val d2 = d*d
val d3 = d2*d
val d4 = d3*d
val d5 = d4*d
val d6 = d5*d
//Compute lat & lon convert to degree
var lat =
phi1 - nu1 * tan(phi1)/rho1 *
(d2/2.0 - (5.0 + 3.0*t1 + 10*c1 - 4.0*c12 -9.0*ecc2Sq) * d4/24.0 +
(61.0 + 90.0*t1 +298.0*c1 + 45.0*t12 -252*ecc2Sq - 3.0*c12) * d6/720.0)
lat = Math.toDegrees(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) / cos(phi1)
lon = Math.toDegrees(lon)
return LatLng(lat,lon)
【讨论】：
以上是关于如何在 python 或 Javascript 中将 UTM 转换为 LatLng的主要内容，如果未能解决你的问题，请参考以下文章