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

Question

我想写一个程序，把纽约市的x/y坐标转换成lat/lng小数点。我是平面/全球映射的新手。我把纽约市在他们网站上提供的常数也包括在内。如果有一篇关于如何做到这一点的好文章，我很想学习！下面是我写的程序和底部的注释输出，以及理想值应该是什么。我只是在这件事上在黑暗中蹒跚而行。

#!/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

我被卡住了，我有很多需要地理编码的东西 谢谢你的帮助！

Answer 1

这些公式可以帮助您：

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

Answer 2

一个字的答案：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])

Answer 3

噢。你最好用图书馆来做这个。稍微搜索一下，应该是the python interface to gdal

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

你最好在gis stackexchange询问更多信息。

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