Конвертировать плоские х / у в лат / лонг - PullRequest
Новая
       60

Конвертировать плоские х / у в лат / лонг

2 голосов
/ 27 марта 2012

Я пытаюсь написать программу, которая принимает координаты x / y города Нью-Йорка и превращает их в десятичные точки широты и долготы.Я новичок в планарном / глобальном отображении.Я включил константы, которые 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

Я застрял, у меня есть тонна таких, которые нуждаются в геокодировании. Спасибо за любую помощь!

Ответы [ 4 ]

5 голосов
/ 28 марта 2012

Ответ одним словом: 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])
1 голос
/ 28 марта 2012
0 голосов
/ 28 марта 2012

Вместо того, чтобы пытаться проработать всю математику, вы можете просто выбрать сетку на поверхности вашей карты и узнать широту / длину этих точек сетки, а затем использовать интерполяцию для выполнения преобразования.В зависимости от линейности проекции для получения хорошей точности может не потребоваться много точек.

0 голосов
/ 28 марта 2012

Owww.вам лучше использовать библиотеку для этого.небольшой поиск показывает, что должен быть интерфейс Python для GDAL

этот вопрос использует GDAL, но не через API Python (они просто вызывают GDAL через командную строкуиз Python), но может помочь.

вам лучше всего спросить gis stackexchange для получения дополнительной информации.

Мне неясно, откуда вы взяли код выше,если вы ссылаетесь на него, я / кто-то может проверить очевидные ошибки реализации.

...