Konwersja współrzędnych X, Y na długość / długość za pomocą pyproj i Proj.4 zwraca nieprawidłowe współrzędne

10

Piszę skrypt w języku Python, który odczytuje wiele plików XML zawierających współrzędne xiy i łączy je wszystkie w jeden plik csv. Szerokość i długość geograficzna są polami wymaganymi w pliku csv, ale mam trudności z konwersją współrzędnych x, y w Ohio North State Plane usFt na 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)

Obie powyższe metody zwracają ten sam wynik, jednak ten długi czas jest gdzieś w Zatoce Hudsona. Kiedy wykreślam współrzędne w ArcMap, prawidłowa długość to: -81.142311,41.688205.

* Zwróć uwagę, że wszystkie długości są dostarczane długo, ponieważ jest to kolejność stosowana przez Proj

Czy ktoś wie, dlaczego otrzymuję nieprawidłowe współrzędne z Proj.4 i pyproj?

Brian
źródło

Odpowiedzi:

8

Otrzymuję takie same wyniki jak @geographika, gdy uruchamiam gdaltransformi narzędzie 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

Odwrócenie współrzędnych xiy punktu daje jednak wynik widoczny w ArcMap:

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

Musisz więc sprawdzić wzrokowo, aby upewnić się, że współrzędne xiy są odpowiednio odwrócone. To był problem, który miałem w przeszłości, kiedy otrzymujesz dwa wyniki, które są wystarczająco podobne, sprowadzasz to do błędu zaokrąglania lub coś takiego.

MerseyViking
źródło
19

PyProj zakłada, że ​​twoje współrzędne są w metrach. Myślę, że przyczyną problemu jest coś związanego ze stopami / metrami.

Wywołanie instancji klasy Proj z argumentami lon, lat przekonwertuje lon / lat (w stopniach) na x / y natywne współrzędne rzutowania mapy (w metrach)

Jeśli opcjonalne słowo kluczowe „preserve_units” ma wartość True, jednostki we współrzędnych rzutowania mapy nie są zmuszane do liczenia.

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

Czy twoje początkowe współrzędne w stopach? Kiedy ładujesz dane do ArcMap, jakich jednostek używa mapa?

To przybliża współrzędne:

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)

Podobny problem można znaleźć tutaj .

geografia
źródło
Dziękuję Ci bardzo. Argument preserve_units zdecydowanie zadziałał, ale współrzędne są nadal niepoprawne. @MerseyViking Ta odpowiedź dała mi prawidłowe współrzędne. Chciałbym oznaczyć obie odpowiedzi jako odpowiedź, ponieważ obie pomogły.
Brian
Cóż, jeśli ludzie głosują na odpowiedź @ geographika bardziej niż moja, wszystko się rozwinie :) Cieszę się, że wszystko działało.
MerseyViking
ponieważ link jest zepsuty, pomocne może być pokazanie, że możesz napisać:p1 = Proj( init="epsg:3734", preserve_units=True )
BenjaminGolder
4

W rzeczywistości próbowałem zrobić to samo, z wyjątkiem siatki stanu południowego OH i natknąłem się na twoje pytanie. Miałem złe wyniki z 3735, teraz otrzymuję prawidłowe wyniki z 3729. Spodziewam się, że jeśli zmienisz z 3734 na 3728, otrzymasz poprawne wyniki.

EPSG: 3728: NAD83 (NSRS2007) / Ohio North (ftUS) EPSG: 3729: NAD83 (NSRS2007) / Ohio South (ftUS) EPSG: 3734: NAD83 / Ohio North (ftUS) EPSG: 3735: NAD83 / Ohio South (ftUS)

Użyłem twojego zapewnionego lat, długi i jestem mniej niż o jedną stopę.

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

p2 (-81.142311,41.688205)

(2339326.6558868014, 739401.4226131936)

vs 2339327,3, 739400,91

inżynier górnictwa
źródło