Das Konvertieren von X-, Y-Koordinaten in Lat / Long mit pyproj und Proj.4 gibt die falschen Koordinaten zurück

10

Ich schreibe ein Python-Skript, das mehrere XML-Dateien mit x- und y-Koordinaten liest und sie alle zu einer einzigen CSV-Datei kombiniert. Breiten- und Längengrade sind Pflichtfelder in der CSV, aber ich habe Schwierigkeiten, die x, y-Koordinaten in der Ohio North State Plane usFt in WGS84 umzuwandeln.

>>> 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)

Beide oben genannten Methoden liefern das gleiche Ergebnis, jedoch befindet sich dieser Lat Long irgendwo in Hudson Bay. Wenn ich die Koordinaten in ArcMap zeichne, lautet die korrekte Lat-Länge: -81,142311,41,688205.

* Beachten Sie, dass alle Lat Longs Long bereitgestellt werden, da dies die Reihenfolge ist, die Proj verwendet

Weiß jemand, warum ich von Proj.4 und pyproj die falschen Koordinaten bekommen würde?

Brian
quelle

Antworten:

8

Ich erhalte die gleichen Ergebnisse wie @geographika, wenn ich gdaltransformund das proj.4-Tool ausführen 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

Das Umkehren der x- und y-Koordinaten Ihres Punkts führt jedoch zu dem Ergebnis, das Sie in ArcMap sehen:

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

Sie müssen also eine visuelle Überprüfung durchführen, um sicherzustellen, dass Ihre x- und y-Koordinaten richtig herum sind. Es ist ein Problem, das ich in der Vergangenheit hatte, wenn Sie zwei Ergebnisse erhalten, die ähnlich genug sind, dass Sie es auf Rundungsfehler oder ähnliches zurückführen.

MerseyViking
quelle
19

PyProj geht davon aus, dass Ihre Koordinaten in Metern angegeben sind. Ich würde vermuten, dass etwas in Bezug auf Fuß / Meter die Ursache des Problems ist.

Wenn Sie eine Proj-Klasseninstanz mit den Argumenten lon aufrufen, konvertiert lat lon / lat (in Grad) in native x / y-Kartenprojektionskoordinaten (in Metern).

Wenn das optionale Schlüsselwort 'keep_units' True ist, müssen die Einheiten in Kartenprojektionskoordinaten keine Meter sein.

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

Sind Ihre Anfangskoordinaten in Fuß? Welche Einheiten verwendet die Karte, wenn Sie die Daten in ArcMap laden?

Dies bringt die Koordinaten etwas näher:

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)

Ein ähnliches Problem finden Sie hier .

geographika
quelle
Vielen Dank. Das Argumenterve_units hat definitiv den Trick gemacht, aber die Koordinaten sind immer noch falsch. @MerseyViking Diese Antwort gab mir die richtigen Koordinaten. Ich wünschte, ich könnte beide Antworten als Antwort markieren, weil beide geholfen haben.
Brian
Nun, wenn die Leute die Antwort von @ geographika mehr als meine positiv bewerten, wird alles gut gehen :) Ich bin froh, dass alles funktioniert hat.
MerseyViking
Da der Link defekt ist, kann es hilfreich sein zu zeigen, dass Sie schreiben können:p1 = Proj( init="epsg:3734", preserve_units=True )
BenjaminGolder
4

Ich habe tatsächlich versucht, dasselbe zu tun, außer mit dem OH South State Plane Grid, und bin auf Ihre Frage gestoßen. Ich habe mit 3735 falsche Ergebnisse erzielt, jetzt erhalte ich mit 3729 die richtigen Ergebnisse. Ich gehe davon aus, dass Sie die richtigen Ergebnisse erhalten, wenn Sie von 3734 auf 3728 wechseln.

EPSG: 3728: NAD83 (NSRS2007) / Ohio Nord (ftUS) EPSG: 3729: NAD83 (NSRS2007) / Ohio Süd (ftUS) EPSG: 3734: NAD83 / Ohio Nord (ftUS) EPSG: 3735: NAD83 / Ohio Süd (ftUS)

Ich habe Ihren bereitgestellten Latz verwendet, bin lang und bin weniger als einen Fuß entfernt.

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

p2 (-81,142311,41,688205)

(2339326.6558868014, 739401.4226131936)

vs 2339327.3, 739400.91

Bergbauingenieur
quelle