Perbedaan lokasi tujuan antara pyproj dan geopy

8

Saya mencari cara untuk menghitung (dengan Python) lokasi tujuan dari titik yang diberikan bantalan dan jangkauan.

Berdasarkan hasil perbandingan dari 2 perpustakaan dalam subjek ( geopydan pyproj), saya perhatikan bahwa ada perbedaan yang meningkat dalam hasil akhir. Sebagai contoh, pada 100 km kira-kira dari urutan desimeter. Ini adalah contoh minimal dari apa yang saya maksud:

from __future__ import absolute_import, division, print_function

long_1 = -1.729722
lat_1 = 53.320556
bearing = 96.021667
distance = 124.8  #in km

# using geopy

import geopy
from geopy.distance import VincentyDistance

origin = geopy.Point(lat_1, long_1)
destination = VincentyDistance(kilometers=distance).destination(origin, bearing)
gp_lat_2 = destination.latitude
gp_long_2 = destination.longitude

# using pyproj

from pyproj import Geod
g = Geod(ellps='WGS84')
prj_long_2, prj_lat_2, prj_bearing_2 = g.fwd(long_1, lat_1, bearing, distance*1000)

# results comparison
print("        | pyproj        | geopy")
print("long_2   %.6f     %.6f" % (prj_long_2, gp_long_2))
print("lat_2    %.6f     %.6f" % (prj_lat_2, gp_lat_2))

print("> DELTA pyproj, geopy")
print("delta long_2   %.7f" % (prj_long_2 - gp_long_2))
print("delta lat_2    %.7f" % (prj_lat_2 - gp_lat_2))

Saya mendapat hasil ini:

        | pyproj        | geopy
long_2   0.127201         0.127199
lat_2    53.188432        53.188432

> DELTA pyproj, geopy
delta long_2   0.0000021
delta lat_2    -0.0000002

Pertanyaan utama saya adalah apakah saya melakukan sesuatu yang salah (kedua pengaturan harus menggunakan WGS84).

Jika tidak, apakah perbedaannya karena perbedaan formula yang digunakan (Vincenty for geopyvs Karney for pyproj)? Misalnya, kesalahan bulat dikutip di sini .

Gmas80
sumber

Jawaban:

6

Sepertinya Anda sudah melakukan semuanya dengan benar. Anda dapat mengevaluasi kesalahan dari setiap metode dengan melakukan perhitungan terbalik untuk menemukan jarak yang diberikan koordinat asal dan tujuan, kemudian mengevaluasi residu jarak. Ini adalah latihan bolak-balik.

# For Vincenty's method:
geopy_inv_dist = geopy.distance.vincenty(origin, destination).m
# For Karney's method:
prj_inv_dist = g.inv(long_1, lat_1, prj_long_2, prj_lat_2)[2]  # s12

print("> inverse distance residule (m)")
print("geopy  %.7f" % (distance * 1000 - geopy_inv_dist))
print("prj    %.7f" % (distance * 1000 - prj_inv_dist))

Menunjukkan:

> inverse distance residule (m)
geopy  0.1434377
prj    0.0000000

Jadi, Anda dapat melihat bahwa metode Vincenty menentukan jarak terbalik yang melebihi desimal yang berbeda untuk koordinat yang sama. Metode Karney memiliki kesalahan dalam ketelitian alat berat, yaitu kurang dari 15 nanometer. Dalam contoh ini kesalahannya adalah 0,1455 nm, yaitu sekitar diameter atom hidrogen.


Masalahnya mungkin dengan metode tujuan geopy . Mari kita bandingkan implementasi kedua metode Vincenty dengan PostGIS versi 2.1, yang ditunjukkan di sini . (PostGIS versi 2.2 dengan Proj 4.9 dan kemudian menggunakan metode Karney ). Residu jarak pulang pergi dari PostGIS 2.1 selalu kurang dari 1 cm. Untuk contoh ini adalah 255 nm:

SELECT PostGIS_Version(),
  ST_AsText(origin) AS origin,
  ST_AsText(ST_Project(origin, distance, azimuth)) AS destination,
  ST_Distance(ST_Project(origin, distance, azimuth), origin) AS roundtrip_distance,
  distance - ST_Distance(ST_Project(origin, distance, azimuth), origin) AS postgis_residual
FROM (
  SELECT 124.8 * 1000 AS distance, radians(96.021667) AS azimuth,
    ST_SetSRID(ST_MakePoint(-1.729722, 53.320556), 4326)::geography AS origin
) AS f;
-[ RECORD 1 ]------+-----------------------------------------
postgis_version    | 2.1 USE_GEOS=1 USE_PROJ=1 USE_STATS=1
origin             | POINT(-1.729722 53.320556)
destination        | POINT(0.12720134063267 53.1884316458524)
roundtrip_distance | 124799.999999745
postgis_residual   | 2.54993210546672e-007
Mike T
sumber
Terima kasih atas jawaban yang sangat teliti dan sangat cerdas ini!
Jason Scheirer
@ mike-t, terima kasih atas jawaban yang bagus! Saya menyukai ide sumber ketiga .. Saya telah membuka geopytiket: github.com/geopy/geopy/issues/140
gmas80