평면 x/y를 위도/경도로 변환

Question

뉴욕시 x/y 좌표를 가져 와서 위/아래 소수점으로 변환하는 프로그램을 작성하려고합니다. 새로운 평면/글로브 매핑. 필자는 NYC가 웹 사이트에 제공 한 상수를 포함 시켰습니다. 또한이 작업을 수행하는 방법에 대한 좋은 기사가 있으면 배우고 싶습니다! 아래는 내가 밑줄에 주석 처리 된 결과와 함께 작성한 프로그램과 이상적인 값이 있어야하는 프로그램입니다. Im는 이것 위에서 어두운 데에서 단지 비틀 거리는 단지 약간을 만난다.

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

한 마디 대답 : 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 2

owww. 이 라이브러리를 사용하는 것이 더 낫습니다. 약간의 검색은 그 the python interface to gdal

this question이 GDAL를 사용하지만 파이썬 API를 통해 (그들은 단지 파이썬 내에서 명령 줄을 통해 GDAL를) 호출하지해야하지만, 도움이 될 수 있습니다 제안합니다.

자세한 내용은 gis stackexchange에 문의하는 것이 가장 좋습니다.

위 코드를 어디서 얻었는지 확실하지 않습니다. 만약 당신이 그것에 링크를 나는/누군가가 명백한 구현 오류를 확인할 수 있습니다.

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

Answer 3

이 공식은 당신을 도움이 될 것입니다 그런 다음 보간법을 사용하여 변환을 수행하십시오. 투영의 선형성에 따라 좋은 정확도를 얻으려면 많은 포인트를 차지하지 않을 수 있습니다.

Answer 4

는 오히려 모든 수학을 통해 해결하려고 노력하는 것보다, 당신은 단지지도의 표면에 격자를 선택할 수 있고 위도/경도 그 그리드의 발견 :