Преобразование координат X, Y в широту / долготу с использованием pyproj и Proj.4 возвращает неправильные координаты

10

Я пишу скрипт на python, который читает несколько файлов XML, содержащих координаты x и y, и объединяет их все в один файл CSV. Широта и долгота являются обязательными полями в CSV, но у меня возникли трудности с преобразованием координат x, y в плоскости штата Северный штат Огайо usFt в WGS84.

>>> p = Proj(r'+proj=lcc +lat_1=41.7 +lat_2=40.43333333333333 +lat_0=39.66666666666666 +lon_0=-82.5 +x_0=600000 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=us-ft +no_defs') #Nad83 State Plane Ohio North US Feet Proj object using parameters
>>> p(739400.91,2339327.3,inverse=True)
(-80.138057868777224, 60.278230707978487)
>>> p1 = Proj(init="epsg:3734") #Nad83 State Plane Ohio North US Feet Proj object using EPSG code
>>> p1(739400.91,2339327.3,inverse=True)
(-80.138057868777224, 60.278230707978487)

Оба метода, приведенные выше, возвращают один и тот же результат, но этот последний длинный где-то в Гудзоновом заливе. Когда я строю координаты в ArcMap, правильный лат длинный: -81.142311,41.688205.

* Обратите внимание, что все значения long указаны в long, lat, так как это порядок, используемый Proj

Кто-нибудь знает, почему я получаю неправильные координаты от Proj.4 и pyproj?

Брайан
источник

Ответы:

8

Я получаю те же результаты, что и @geographika, когда запускаю gdaltransformи инструмент proj.4 cs2cs:

$ gdaltransform -s_srs EPSG:3734 -t_srs EPSG:4326
739400.9 2339327.3             
-87.3195485720169 45.9860670658218 0

cs2cs +proj=lcc +lat_1=41.7 +lat_2=40.43333333333333 +lat_0=39.66666666666666 +lon_0=-82.5 +x_0=600000 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=us-ft +no_defs +to +proj=lonlat +datum=WGS84
739400.9 2339327.3
87d19'10.375"W 45d59'9.841"N 0.000

Однако изменение координат x и y вашей точки дает результат, который вы видите в ArcMap:

gdaltransform -s_srs EPSG:3734 -t_srs EPSG:4326
2339327.3 739400.9
-81.1423086719059 41.6882035384526 0

Так что вам нужно будет сделать визуальную проверку, чтобы убедиться, что ваши координаты x и y правильные и правильные. Это проблема, с которой я сталкивался в прошлом, когда вы получаете два результата, которые достаточно похожи, вы объясняете это ошибкой округления или чем-то еще.

MerseyViking
источник
19

PyProj предполагает, что ваши координаты в метрах. Я думаю, что что-то, связанное с футами / метрами, является причиной проблемы.

Вызов экземпляра класса Proj с аргументами lon, lat преобразует lon / lat (в градусах) в координаты проекции собственной карты x / y (в метрах)

Если для необязательного ключевого слова preserve_units задано значение True, единицы измерения в координатах проекции карты не обязательно должны быть метрами.

http://pyproj.googlecode.com/svn/trunk/docs/pyproj.Proj-class.html

Ваши начальные координаты в футах? Когда вы загружаете данные в ArcMap, какие единицы измерения использует карта?

Это немного приближает координаты:

p1 = Proj(init="epsg:3734")
#1 foot = 0.3048 meters
conv = 0.3048
print p1(739400.91 * conv,2339327.3 * conv,inverse=True)
(-87.3195533069909, 45.98605408134072)

Подобную проблему можно найти здесь .

geographika
источник
Большое спасибо. Аргумент preserve_units определенно добился цели, но координаты все еще неверны. @MerseyViking Этот ответ дал мне правильные координаты. Я хотел бы отметить оба ответа как ответ, потому что они оба помогли.
Брайан
Ну, если люди ответят больше, чем мой ответ, то все получится :) Рад, что все сработало.
MerseyViking
поскольку ссылка не работает, было бы полезно показать, что вы можете написать:p1 = Proj( init="epsg:3734", preserve_units=True )
BenjaminGolder
4

На самом деле я пытался сделать то же самое, за исключением сетки самолетов южного штата ОХ, и я столкнулся с вашим вопросом. Я получал неправильные результаты с 3735, теперь я получаю правильные результаты с 3729. Я ожидаю, что если вы измените с 3734 на 3728, вы получите правильные результаты.

EPSG: 3728: NAD83 (NSRS2007) / Север Огайо (ftUS) EPSG: 3729: NAD83 (NSRS2007) / Юг Огайо (ftUS) EPSG: 3734: NAD83 / Север Огайо (ftUS) EPSG: 3735: NAD83 / Юг Огайо (ftUS)

Я использовал предоставленный вам широчайший, длинный и ушел менее чем на один фут.

p2 = pyproj.Proj (init = "epsg: 3728", preserve_units = True)

р2 (-81.142311,41.688205)

(2339326.6558868014, 739401.4226131936)

против 2339327,3, 739400,91

горный инженер
источник