Was ist die Einheit des Formlängenattributs?

11

Ich mache eine sehr einfache Berechnung der Länge einer Polylinie mit Shapely:

from shapely.geometry import LineString
... 
xy_list = [map(float,e) for e in xy_intm]
line = LineString(xy_list)
s = '%s,%s,%s' % (fr,to,line.length)

Meine Koordinaten sind in WGS84. Ich kann anscheinend keine Informationen über das Längenattribut von Shapely finden. Was ist die Einheit des Längenattributs? Gibt es eine einfache Möglichkeit, in km oder Meter umzurechnen?

LarsVegas
quelle
Können Sie die Koordinaten und die Länge für zwei Beispielformen angeben?
Vince

Antworten:

13

Wie Alfaciano formschön sagt, ist der Abstand der euklidische Abstand oder der lineare Abstand zwischen zwei Punkten auf einer Ebene und nicht der Großkreisabstand zwischen zwei Punkten auf einer Kugel.

from shapely.geometry import Point
import math


point1 = Point(50.67,4.62)
point2 = Point(51.67, 4.64)

# Euclidean Distance
def Euclidean_distance(point1,point2):
     return math.sqrt((point2.x()-point1.x())**2 + (point2.y()-point1.y())**2)

print Euclidean_distance(point1,point2)
1.00019998 # distance in degrees (coordinates of the points in degrees)

# with Shapely
print point1.distance(point2)
1.0001999800039989 #distance in degrees (coordinates of the points in degrees)

Für die Großkreisentfernung, müssen Sie Algorithmen als Kosinussatz oder der Haversine Formel verwenden (siehe Warum ist Kosinussatz stärker bevorzugt als Haversine , wenn der Abstand zwischen zwei Breiten-Längen - Punkte Berechnung? ) Oder das Modul verwenden pyproj dass führt geodätische Berechnungen durch.

# law of cosines
distance = math.acos(math.sin(math.radians(point1.y))*math.sin(math.radians(point2.y))+math.cos(math.radians(point1.y))*math.cos(math.radians(point2.y))*math.cos(math.radians(point2.x)-math.radians(point1.x)))*6371
print "{0:8.4f}".format(distance)
110.8544 # in km
# Haversine formula
dLat = math.radians(point2.y) - math.radians(point1.y)
dLon = math.radians(point2.x) - math.radians(point1.x)
a = math.sin(dLat/2) * math.sin(dLat/2) + math.cos(math.radians(point1.y)) * math.cos(math.radians(point2.y)) * math.sin(dLon/2) * math.sin(dLon/2)
distance = 6371 * 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
print "{0:8.4f}".format(distance)distance
110.8544 #in km

# with pyproj
import pyproj
geod = pyproj.Geod(ellps='WGS84')
angle1,angle2,distance = geod.inv(point1.x, point1.y, point2.x, point2.y)
print "{0:8.4f}".format(distance/1000)
110.9807 #in km

Sie können das Ergebnis mit dem Longitude Latitude Distance Calculator testen

Geben Sie hier die Bildbeschreibung ein

Gen
quelle
Tolle Antwort, Gene! Vielen Dank für Ihre sehr ausführliche Erklärung.
Antonio Falciano
1
In der Tat eine gute Antwort. Wenn ich mich nicht irre, gibt es ein anderes Python-Paket namens geopy, das eine Großkreisentfernung und eine Vincenty-Entfernungsberechnung implementiert hat.
LarsVegas
Hier einige Details zur Berechnung der geodätischen Entfernung mit geopy.
Antonio Falciano
13

Koordinatensystem

[...] Shapely unterstützt keine Koordinatensystemtransformationen. Alle Operationen an zwei oder mehr Merkmalen setzen voraus, dass die Merkmale in derselben kartesischen Ebene vorhanden sind.

Quelle: http://toblerity.org/shapely/manual.html#coordinate-systems

Da es shapelyin Bezug auf SRS völlig agnostisch ist, ist es ziemlich offensichtlich, dass das Längenattribut in derselben Koordinateneinheit Ihres Linienstreifens ausgedrückt wird, dh in Grad. Eigentlich:

>>> from shapely.geometry import LineString
>>> line = LineString([(0, 0), (1, 1)])
>>> line.length
1.4142135623730951

Wenn Sie stattdessen die Länge in Metern ausdrücken möchten, müssen Sie Ihre Geometrien mithilfe von pyproj von WGS84 in einen projizierten SRS umwandeln (oder besser eine geodätische Entfernungsberechnung durchführen, siehe Gens Antwort). Im Detail werden seit Version 1.2.18 ( shapely.__version__) shapelydie Geometrie-Transformationsfunktionen ( http://toblerity.org/shapely/shapely.html#module-shapely.ops ) unterstützt, mit denen wir sie verwenden können pyproj. Hier ist ein kurzes Beispiel:

from shapely.geometry import LineString
from shapely.ops import transform
from functools import partial
import pyproj

line1 = LineString([(15.799406, 40.636069), (15.810173,40.640246)])
print str(line1.length) + " degrees"
# 0.0115488362184 degrees

# Geometry transform function based on pyproj.transform
project = partial(
    pyproj.transform,
    pyproj.Proj(init='EPSG:4326'),
    pyproj.Proj(init='EPSG:32633'))

line2 = transform(project, line1)
print str(line2.length) + " meters"
# 1021.77585965 meters
Antonio Falciano
quelle