经纬度转UTM再转回经纬度极为错误，为什么？

Question

我尝试过以下方法，

输入的是经纬度数据，

然后我会计算一个围绕它的框，假设是50米，所以在东西和南北方向上各加减50米。

现在我再把它转换回经纬度，并用一个脚本：

http://robotics.ai.uiuc.edu/~hyoon24/LatLongUTMconversion.py，结果却完全不对，之前的经度大约是7，之后变成了大约2。

zone, easting, northing = LLtoUTM(23, location.get_lat(), location.get_lon()) 

topUTM = northing + error
bottomUTM = northing - error
leftUTM = easting - error
rightUTM = easting + error
left, top = UTMtoLL(23, leftUTM, topUTM, zone)

是我的代码有问题，还是这个脚本本身有缺陷呢？

所以我尝试使用pyproj，只是把经纬度转换成utm再转换回经纬度，看看会发生什么。

>>> p = pyproj.Proj(proj='utm', zone=32, ellps='WGS84')
>>> p
<pyproj.Proj object at 0x7ff9b8487dd0>
>>> x,y = p(47.9941214, 7.8509671)
>>> print x,y
5159550.36822 1114087.43925
>>> print p(x,y,inverse=True)
(47.971558538495991, 7.8546573140162605)

结果没有上面的脚本偏差那么大，但还是有明显的错误，无法使用。这是为什么呢？我该怎么做才能得到更准确的结果？

编辑：

我运行了test()，所有测试都通过了。

在epsg文件中没有这样的东西。我找到的最接近的是这个：

<32632> +proj=utm +zone=32 +ellps=WGS84 +datum=WGS84 +units=m +no_defs <>

没有tmerc。此外，我需要传递towgs84作为参数吗？就是上面那些吗？

Answer 1

我对pyproj没有任何问题，可以试试下面的代码

from pyproj import Proj

Lat = 52.063098675
Lon = -114.132980348 #Calgary

ZoneNo = "11" #Manually input, or calcuated from Lat Lon
myProj = Proj("+proj=utm +zone="+ZoneNo+",\
+north +ellps=WGS84 +datum=WGS84 +units=m +no_defs") #north for north hemisphere
UTMx, UTMy = myProj(Lon, Lat)

########################################

#UTM ==> Lat Lon:
ZoneNo = "11" #Manually input or from other sources
myProj = Proj("+proj=utm +zone="+\
ZoneNo+", +north +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
Lon2, Lat2 = myProj(UTMx, UTMy,inverse=True)

print Lat2
print Lon2

Answer 2

你的代码有问题。

首先，之前提到的PyProj问题确实存在。你应该检查一下你的epsg文件，确保里面有这一行：

<2392> +proj=tmerc +lat_0=0 +lon_0=24 +k=1.000000 +x_0=2500000 +y_0=0 +ellps=intl +towgs84=-90.7,-106.1,-119.2,4.09,0.218,-1.05,1.37 +units=m +no_defs no_defs <>

注意一下 towgs84 这个参数。

你在使用PyProj时出错是因为命令用错了。

比如说，如果我们把47.9941214N, 7.8509671E这个坐标转换成UTM坐标，结果是区域32，414278东坐标，5316286北坐标。

你执行了以下的PyProj操作：

p = pyproj.Proj(proj='utm', zone=32, ellps='WGS84')
>>> x,y = p(47.9941214, 7.8509671)
>>> print x,y
5159550.36822 1114087.43925
>>> print p(x,y,inverse=True)
(47.971558538495991, 7.8546573140162605)

但是，如果我们查看PyProj的文档，会看到以下内容：

用经度和纬度调用一个Proj类实例，会把经度/纬度（以度为单位）转换成x/y本地地图投影坐标（以米为单位）。

我们再试着运行一下你的PyProj操作，不过这次把经度和纬度的顺序调换一下：

p = pyproj.Proj(proj='utm', zone=32, ellps='WGS84')
>>> x,y = p(7.8509671, 47.9941214)
>>> print x,y
414278.16731 5316285.59492
>>> print p(x,y,inverse=True)
(7.850967099999812, 47.994121399999784)

这样操作几乎完美地反转了结果！

关于你问题的第一部分，如果你查看一下http://robotics.ai.uiuc.edu/~hyoon24/LatLongUTMconversion.py中UTMtoLL的定义，你会发现：

UTMtoLL(ReferenceEllipsoid, northing, easting, zone)

然而，你却用UTMtoLL(23, leftUTM, topUTM, zone)，其中leftUTM是东坐标，topUTM是北坐标。

所以在你第一个脚本和PyProj的情况下，你都用了错误的参数顺序。

这提醒我们在指责别人之前，最好先仔细检查自己的工作。话虽如此，Python的文档并不是特别好，而且在这个情况下，PyProj的文档也很难懂。如果有一个更清晰的网页解释和使用示例，可能就能避免你现在的困惑了。

Answer 3

我上周创建了一个小的UTM转换库，用于Python，并把它上传到了Python包索引网站：http://pypi.python.org/pypi/utm

我把它和pyproj进行了比较，发现这个库更快，准确度也更高。根据你提供的样本数据，结果如下：

>>> import utm

>>> u = utm.from_latlon(47.9941214, 7.8509671)
>>> print u
(414278, 5316285, 32, 'T')

>>> print utm.to_latlon(*u)
(47.994157948891505, 7.850963967574302)

更新：Richards在下面的回答中描述了这个问题的真正解决方案。