经纬度转UTM再转回经纬度的方法存在极大缺陷，为什么？

Question

经纬度转UTM再转回经纬度的方法存在极大缺陷，为什么？

pythonlatitude-longitudeutm

19

我尝试了以下步骤：

输入：经纬度数据

然后，我会计算一个以该点为中心，东/北坐标值分别加减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？就上面那些参数吗？

- luh

4个回答

34

错误在您的代码中。

首先，另一个答案中列出的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文档，可以看到如下内容:

使用经度和纬度(lon, lat)调用一个Proj类实例将把经度/纬度(以度为单位)转换为x/y本地地图投影坐标(以米为单位)。

让我们再次尝试运行OP的PyProj操作，但是交换lon/lat参数的顺序:

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的文档最多也只是晦涩难懂。一个漂亮的基于Web的解释和伴随使用示例的命令可能已经防止了你的烦恼。

- Richard

1

有没有办法从经纬度自动检测UTM带号？ - BetterCallMe

自2013年以来，这方面有更好的文档吗？（GDAL似乎在往返经纬度-> UTM-> 经纬度时出现问题，请参见此[图表]（https://github.com/corteva/rioxarray/issues/396）） - denis

@Lostman：math.floor(((lon + 180) % 360) / 6) + 1应该可以。请参见https://dev59.com/DGox5IYBdhLWcg3wekIr - Eric Duminil

4

我对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

- Ting On Chan

2

您在pyProj上遇到的问题，听起来很像这里描述的问题：http://code.google.com/p/pyproj/issues/detail?id=3 ，已经被解决了： solved! in epsg file there must be

<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，请查看该线程。此外，模块的test()函数是否有效？您是否尝试过带有test目录中的脚本？

- agf

那些希望Ubuntu 13.04 Raring包含最新的PyProj的人会发现自己是错误的。需要从源代码安装！ - Richard

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- TBieniek · Accepted Answer

上周我创建了一个小型的Python UTM转换库，并将其上传到Python Package Index：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在下面的答案中描述了这个问题的真正解决方案。