Qu'est-ce que l'attribut unité de longueur galbée?

11

Je fais un calcul très simple de la longueur d'une polyligne en utilisant galbée:

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)

Mes coordonnées sont en WGS84. Je n'arrive pas à trouver d'informations sur l'attribut de longueur galbée. Quelle est l'unité de l'attribut de longueur? Existe-t-il un moyen facile de convertir en km ou en mètres?

LarsVegas
la source
Pouvez-vous fournir les coordonnées et la longueur de deux exemples de formes?
Vince

Réponses:

13

Comme le dit alfaciano dans galbée , la distance est la distance euclidienne ou la distance linéaire entre deux points sur un plan et non la distance du grand cercle entre deux points sur une sphère.

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)

Pour la distance du grand cercle, vous devez utiliser des algorithmes comme loi des cosinus ou la formule Haversine (regardez pourquoi la loi des cosinus est-elle plus préférable que la haversine lors du calcul de la distance entre deux points de latitude-longitude? ) Ou utilisez le module pyproj qui effectue des calculs géodésiques.

# 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

Vous pouvez tester le résultat sur le Calculateur de distance de Latitude Latitude

entrez la description de l'image ici

gène
la source
Excellente réponse, Gene! Merci beaucoup pour votre explication très détaillée.
Antonio Falciano
1
En effet, bonne réponse. Si je ne me trompe pas, il existe un autre package python appelé geopyqui a implémenté un calcul de distance grand cercle et de distance Vincenty.
LarsVegas
Voici quelques détails sur le calcul de la distance géodésique avec geopy.
Antonio Falciano
13

Systèmes de coordonnées

[...] Shapely ne prend pas en charge les transformations du système de coordonnées. Toutes les opérations sur deux ou plusieurs entités supposent que les entités existent dans le même plan cartésien.

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

Étant shapelycomplètement agnostique par rapport à SRS, il est tout à fait évident que l'attribut de longueur est exprimé dans la même unité de coordonnées de votre chaîne de lignes, à savoir les degrés. En réalité:

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

Au lieu de cela, si vous souhaitez exprimer la longueur en mètres, vous devez transformer vos géométries du WGS84 en un SRS projeté en utilisant pyproj (ou, mieux, exécuter le calcul de la distance géodésique, voir la réponse de Gene). En détail, depuis la version 1.2.18 ( shapely.__version__), shapelyprend en charge les fonctions de transformation de la géométrie ( http://toblerity.org/shapely/shapely.html#module-shapely.ops ) avec lesquelles nous pouvons l'utiliser en conjonction pyproj. Voici un petit exemple:

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
la source