我正在尝试编写一个程序，将纽约市的 x/y 坐标转换为纬度/经度小数点。我是平面/全局映射的新手。我包括纽约市在其网站上提供的常量。另外，如果有一篇关于如何做到这一点的好文章，我很乐意学习!下面是我编写的程序以及底部的注释输出以及理想值应该是什么。我有点在这方面摸不着头脑。
#!/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
我很困惑，我有很多需要地理编码的东西 谢谢你的帮助!
最佳答案
一句话回答: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])
关于python - 将平面 x/y 转换为纬度/经度，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/9895873/