将平面 x/y 转换为 lat/long答案

【问题标题】：Convert planar x/y to lat/long将平面 x/y 转换为 lat/long
【发布时间】：2012-03-27 18:58:59
【问题描述】：

我正在尝试编写一个程序，该程序采用纽约市 x/y 坐标并将它们转换为 lat/lng 小数点。我是平面/地球映射的新手。我已经包含了 NYC 在其网站上提供的常量。另外，如果有一篇关于如何做到这一点的好文章，我很想学习！下面是我编写的程序以及底部的注释输出以及理想值应该是什么。我有点在这方面的黑暗中跌跌撞撞。

#!/usr/bin/python
from math import *

"""
Supplied by NYC
Lambert Conformal Conic:

    Standard Parallel: 40.666667
    Standard Parallel: 41.033333
    Longitude of Central Meridian: -74.000000
    Latitude of Projection Origin: 40.166667
    False Easting: 984250.000000
    False Northing: 0.000000

"""

x = 981106                      #nyc x coord
y = 195544                      #nyc y coord
a = 6378137                     #' major radius of ellipsoid, map units (NAD 83)
e = 0.08181922146               #' eccentricity of ellipsoid (NAD 83)
angRad = pi/180                 #' number of radians in a degree
pi4 = pi/4                      #' Pi / 4

p0 = 40.166667 * angRad        #' latitude of origin
p1 = 40.666667 * angRad        #' latitude of first standard parallel
p2 = 41.033333 * angRad        #' latitude of second standard parallel
m0 = -74.000000 * angRad       #' central meridian
x0 = 984250.000000             #' False easting of central meridian, map units

m1 = cos(p1) / sqrt(1 - ((e ** 2) * sin(p1) ** 2))
m2 = cos(p2) / sqrt(1 - ((e ** 2) * sin(p2) ** 2))
t0 = tan(pi4 - (p0 / 2))
t1 = tan(pi4 - (p1 / 2))
t2 = tan(pi4 - (p2 / 2))
t0 = t0 / (((1 - (e * (sin(p0)))) / (1 + (e * (sin(p0)))))**(e / 2))
t1 = t1 / (((1 - (e * (sin(p1)))) / (1 + (e * (sin(p1)))))**(e / 2))
t2 = t2 / (((1 - (e * (sin(p2)))) / (1 + (e * (sin(p2)))))**(e / 2))
n = log(m1 / m2) / log(t1 / t2)
f = m1 / (n * (t1 ** n))
rho0 = a * f * (t0 ** n)

x = x - x0
pi2 = pi4 * 2
rho = sqrt((x ** 2) + ((rho0 - y) ** 2))
theta = atan(x / (rho0 - y))
t = (rho / (a * f)) ** (1 / n)
lon = (theta / n) + m0
x = x + x0

lat0 = pi2 - (2 * atan(t))

part1 = (1 - (e * sin(lat0))) / (1 + (e * sin(lat0)))
lat1 = pi2 - (2 * atan(t * (part1 ** (e / 2))))
while abs(lat1 - lat0) < 0.000000002:
    lat0 = lat1
    part1 = (1 - (e * sin(lat0))) / (1 + (e * sin(lat0)))
    lat1 = pi2 - (2 * atan(t * (part1 ^ (e / 2))))

lat = lat1 / angRad
lon = lon / angRad

print lat,lon
#output : 41.9266666432 -74.0378981653
#should be 40.703778, -74.011829

我很困惑，我有很多需要地理编码的东西谢谢你的帮助！

【问题讨论】：

也许en.wikipedia.org/wiki/Map_projection 可以帮助您入门？
更多关于所用投影的详细信息：en.wikipedia.org/wiki/Lambert_conformal_conic_projection
您的代码表明这是一个圆锥投影，而不是椭圆...
嗯，很好，甚至没看到，这可能就是为什么我只是脱发了。当数据为圆锥时，我正在计算椭球的偏心率。

标签： python coordinates gis geocoding latitude-longitude

【解决方案1】：

一句话回答：pyproj

>>> from pyproj import Proj
>>> pnyc = Proj(
...     proj='lcc',
...     datum='NAD83',
...     lat_1=40.666667,
...     lat_2=41.033333,
...     lat_0=40.166667,
...     lon_0=-74.0,
...     x_0=984250.0,
...     y_0=0.0)
>>> x = [981106.0]
>>> y = [195544.0]
>>> lon, lat = pnyc(x, y, inverse=True)
>>> lon, lat
([-74.037898165369015], [41.927378144152335])

【讨论】：

这看起来棒极了！但仍然不是我需要的准确性，这只是数据在 x,y 中的限制吗？实际上，这只是输出我已经拥有的结果：S
在玩过这个不应该是经度-74之后
我使用 Google 从我拥有的地址进行地理编码。这应该最终在曼哈顿下城的某个地方结束。这是地图描述的链接nyc.gov/html/dcp/html/bytes/…
没关系，地址是曼哈顿的 99 BROAD STREET，即使更改为 -74，我仍然会在波基普西的北部登陆。除非我从 nyc 得到的这些数据是坏的。
现在我真的明白了。纽约市的数据以英尺为单位。 :(

【解决方案2】：

这些公式应该可以帮助您：

http://www.linz.govt.nz/geodetic/conversion-coordinates/projection-conversions/lambert-conformal-conic#lbl3

【讨论】：

【解决方案3】：

哎呀。你最好为此使用图书馆。稍作搜索表明应该是the python interface to gdal

this question 使用 gdal，但不是通过 python api（他们只是通过 python 中的命令行调用 gdal），但可能会有所帮助。

您最好向gis stackexchange 询问更多信息。

我不清楚你从哪里得到上面的代码。如果你链接到它，我/有人可以检查明显的实现错误。

【讨论】：

大家好，谢谢大家！我破解了在蒙大拿州 GIS 存储库中找到的一些代码。它是用 asp/vb 写的，这里有链接。它与我的代码几乎相同。不过我要去看看这个图书馆。 nris.mt.gov/gis/projection/stategeo_asp.txt 另外，我确实有这些建筑物的地址，这就是我获得正确地理编码的方式。我只是对其中许多使用地理编码服务有 waaaaaaaaaaaay。
@arboc7 发现我的数据是圆锥形的，但这个程序是针对椭圆体的。我认为这就是为什么我得到接近答案但不是正确答案的原因。

【解决方案4】：

您可以在地图表面上选择一个网格并找出这些网格点的纬度/经度，然后使用插值进行转换，而不是尝试完成所有数学运算。根据投影的线性度，可能不需要很多点即可获得良好的精度。

【讨论】：