Messabstand im sphärischen Mercator gegen UTM in Zonen

11

Ich habe Punkte in WGS84 lat / long und möchte "kleine" (weniger als 5 km) Entfernungen zwischen ihnen messen.

Ich kann die Haversine-Formel von http://www.movable-type.co.uk/scripts/latlong.html verwenden und sie funktioniert sehr gut.

Ich möchte jedoch Python Shapely-Bibliotheken verwenden, damit ich mehr Operationen als nur Entfernungen ausführen kann, und weil auf der Skala, mit der ich arbeite, eine flache Erde eine ausreichend gute Annäherung ist. Um die geografischen Koordinaten zuverlässig auf eine kartesische Koordinate zu projizieren, verwende ich Pythons proj4, aber es scheinen größere Fehler zu auftreten, als ich möchte.

Wenn ich die lokale UTM-Zone verwende, erhalte ich Unterschiede zwischen dem Haversine von einigen Metern, was in Ordnung ist. Aber ich möchte nicht die UTM-Zone berechnen müssen (die Punkte könnten weltweit sein), also habe ich es mit "sphärischem Mercator" versucht, aber jetzt liegen die Unterschiede zwischen Haversine und projizierten Entfernungen weit über 100%. Ist das wirklich richtig für sphärischen Mercator? Alles, was ich wirklich will, ist eine praktikable kartesische Projektion für zwei Punkte innerhalb von 5 km voneinander auf der ganzen Welt.

from shapely.geometry import Point
from pyproj import Proj

proj = Proj(proj='utm',zone=27,ellps='WGS84')
#proj = Proj(init="epsg:3785")  # spherical mercator, should work anywhere...

point1_geo = (-21.9309694, 64.1455718)
point2_geo = (-21.9372481, 64.1478206)
point1 = proj(point1_geo[0], point1_geo[1])
point2 = proj(point2_geo[0], point2_geo[1])

point1_cart = Point(point1)
point2_cart = Point(point2)

print "p1-p2 (haversine)", hdistance(point1_geo, point2_geo)
print "p1-p2 (cartesian)", point1_cart.distance(point2_cart)

Zu diesem Zeitpunkt beträgt der Abstand zwischen ihnen 394 m und unter Verwendung der utm-Zone 27 395 m. Aber wenn ich sphärischen Mercator benutze, beträgt die kartesische Entfernung 904 m, was weit entfernt ist.

python coordinate-system distance Karl P.
quelle
Die UTM-Zone lässt sich anhand der Längengrade leicht "trainieren". Wählen Sie ein typisches Längengrad-Lambda, -180 <= Lambda <180, und verwenden Sie es, um die Zonennummer als Int ((180 + Lambda) / 6) +1 zu berechnen. Verwenden Sie das Breitengradzeichen, um zwischen Nord und Süd zu entscheiden. Sie müssen die speziellen Polarzonen in hohen Breiten nicht verwenden. In der Tat können Sie fast in der Nähe einer Stange fast jede UTM-Zone verwenden.
whuber

Antworten:

17

Ja, bei einer globalen Mercator-Projektion treten diese Fehler auf: Sie ist am Äquator genau und die Verzerrung nimmt exponentiell mit dem Breitengrad vom Äquator weg zu. Die Abstandsverzerrung beträgt genau 2 (100%) bei 60 Grad Breite. Bei Ihren Testbreiten (64,14 Grad) berechne ich eine Verzerrung von 2,294, was genau mit dem Verhältnis 904/394 = 2,294 übereinstimmt. (Früher habe ich 2.301 berechnet, aber das basierte auf einer Kugel, nicht auf dem WGS84-Ellipsoid. Die Differenz (von 0,3%) gibt uns einen Eindruck von der Genauigkeit, die Sie durch die Verwendung einer ellipsoidbasierten Projektion gegenüber der kugelbasierten Haversine-Formel erzielen können. )

Es gibt keine globale Projektion, die überall hochgenaue Entfernungen liefert. Dies ist einer der Gründe, warum das UTM-Zonensystem verwendet wird!

Eine Lösung besteht darin, für alle Ihre Berechnungen eine sphärische Geometrie zu verwenden, die Sie jedoch abgelehnt haben (was sinnvoll ist, wenn Sie komplexe Operationen ausführen, die Entscheidung jedoch möglicherweise erneut überprüft werden sollte).

Eine andere Lösung besteht darin , die Projektion an die zu vergleichenden Punkte anzupassen . Sie können beispielsweise sicher einen transversalen Mercator (wie im UTM-System) verwenden, dessen Meridian nahe der Mitte des interessierenden Bereichs liegt. Das Verschieben des Meridians ist ganz einfach: Subtrahieren Sie einfach die Länge des Meridians von allen Längen und verwenden Sie eine einzelne TM-Projektion, die auf dem Nullmeridian zentriert ist (mit einem Skalierungsfaktor von 1 anstelle von 0,9996 des UTM-Systems). Für Ihre Arbeit wird dies tendenziell mehr seingenauer als mit UTM selbst. Es liefert korrekte Winkel (TM ist konform) und ist für Punkte, die nur wenige zehn Kilometer voneinander entfernt sind, bemerkenswert genau: Erwarten Sie eine bessere Genauigkeit als sechsstellig. Tatsächlich würde ich eher kleine Unterschiede zwischen diesen angepassten TM-Abständen und den Haversine-Abständen dem Unterschied zwischen dem Ellipsoid (für die TM-Projektion verwendet) und der Kugel (von Haversine verwendet) zuschreiben, als der Verzerrung in der Projektion.

whuber
quelle
Das klingt ziemlich perfekt. Ich denke, ich muss dann meine eigene Init-Zeichenfolge für proj4 erstellen, anstatt eine der vorhandenen EPSG-Zeichenfolgen verwenden zu können.
Karl P
1
+1 Projektion an die Punkte anpassen. Ich bevorzuge Querplatten-Carrée gegenüber Quer-Mercator, aber gegenüber ausreichend kleinen Bereichen ("großer Maßstab") ergibt fast jede "zentrierte" Projektion in der Nähe des interessierenden Bereichs eine gute Genauigkeit.
David Cary
@ David Interessante Idee. Auf der Kugel liegt die Querplatte Carree (Cassini) nahe an der ungefähren Formel, die ich unter gis.stackexchange.com/posts/2964/edit angegeben habe (was hier eine akzeptable Lösung sein könnte). Die Formeln für TM und TPC sind auf der Kugel ähnlich. Auf einem Ellipsoid ist TPC etwas einfacher. TM wird wahrscheinlich von mehr Software unterstützt.
whuber
@Karl Wenn Sie möchten, können Sie jede der TM-Zonen verwenden. Verschieben Sie einfach alle Längengrade so, dass ein zentraler Punkt in Ihrer Region von Interesse mit dem zentralen Meridian der ausgewählten Zone übereinstimmt. Multiplizieren Sie alle Abstände mit 1 / 0,9996 (und multiplizieren Sie alle Flächen mit dem Quadrat dieses Faktors), ändern Sie keine Winkel oder Peilungen und verschieben Sie - wenn Ihre Berechnungen neue Punkte ergeben - einfach ihre Längen zurück zum ursprünglichen Koordinatensystem .
whuber