Apa yang dimaksud dengan unit atribut length bentuk?

11

Saya sedang melakukan perhitungan yang sangat sederhana dari panjang polyline menggunakan rupawan:

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)

Koordinat saya ada di WGS84. Sepertinya saya tidak dapat menemukan informasi tentang atribut panjang shapely. Apa satuan atribut panjang? Apakah ada cara mudah untuk mengkonversi ke km atau meter?

LarsVegas
sumber
Bisakah Anda memberikan koordinat dan panjang untuk dua bentuk sampel?
Vince

Jawaban:

13

Seperti yang dikatakan alfaciano secara rupawan , jaraknya adalah jarak Euclidean atau jarak Linear antara dua titik di pesawat dan bukan jarak lingkaran-besar antara dua titik di bola.

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)

Untuk jarak lingkaran besar, Anda perlu menggunakan algoritma sebagai hukum cosinus atau rumus Haversine (lihat Mengapa hukum kosinus lebih disukai daripada haversine saat menghitung jarak antara dua titik lintang-bujur? ) Atau menggunakan modul pyproj yang melakukan perhitungan geodetik.

# 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

Anda dapat menguji hasilnya di Longitude Latitude Distance Calculator

masukkan deskripsi gambar di sini

gen
sumber
Jawaban yang bagus, Gene! Terima kasih banyak atas penjelasan Anda yang sangat mendetail.
Antonio Falciano
1
Memang, jawaban yang bagus. Jika saya tidak salah ada paket python lain yang disebut geopyyang telah mengimplementasikan jarak lingkaran besar dan penghitungan jarak Vincenty.
LarsVegas
Berikut ini beberapa detail tentang perhitungan jarak geodesik dengan geopy.
Antonio Falciano
13

Sistem Koordinasi

[...] Shapely tidak mendukung transformasi sistem koordinat. Semua operasi pada dua atau lebih fitur menganggap bahwa fitur ada di pesawat Cartesian yang sama.

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

Menjadi shapelysepenuhnya agnostik sehubungan dengan SRS, cukup jelas bahwa atribut panjang dinyatakan dalam satuan koordinat yang sama dengan linestring Anda, yaitu derajat. Faktanya:

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

Sebaliknya, jika Anda ingin mengekspresikan panjang dalam meter, Anda harus mengubah geometri Anda dari WGS84 ke SRS yang diproyeksikan menggunakan pyproj (atau, lebih baik, jalankan perhitungan jarak geodesik, lihat jawaban Gene). Secara terperinci, sejak versi 1.2.18 ( shapely.__version__), shapelymendukung fungsi transformasi geometri ( http://toblerity.org/shapely/shapely.html#module-shapely.ops ) yang dapat kita gunakan bersamaan pyproj. Ini contoh singkatnya:

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
sumber