Pomiar odległości w kulistym Mercatorze względem strefowego UTM

Mam punkty w WGS84 lat / long i chciałbym zmierzyć „małe” (mniej niż powiedzmy 5 km) odległości między nimi.

Mogę użyć formuły haversine z http://www.movable-type.co.uk/scripts/latlong.html i działa bardzo dobrze.

Chciałbym jednak użyć bibliotek Python Shapely, aby móc wykonywać więcej operacji niż tylko odległość, a ponieważ w skali, z którą pracuję, płaska ziemia jest wystarczająco dobrym przybliżeniem. Aby niezawodnie rzutować współrzędne geograficzne na współrzędne kartezjańskie, używam Pythona proj4, ale wydaje się, że dostaję większe błędy, niż bym chciał.

Jeśli korzystam z lokalnej strefy UTM, dostaję różnice między haversine o kilka metrów, co jest w porządku. Ale nie chcę musieć opracowywać strefy UTM (punkty mogą być na całym świecie), więc spróbowałem z „kulistym Mercatorem”, ale teraz różnice między podłużną i przewidywaną odległością są znacznie większe niż 100%. Czy to naprawdę jest odpowiednie dla sferycznego Mercatora? Wszystko, czego naprawdę chcę, to wykonalna projekcja kartezjańska dla dwóch punktów w odległości 5 km od siebie w dowolnym miejscu na świecie.

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)

W tym momencie odległość między nimi wynosi 394 m, a strefa utm 27, 395 m. Ale jeśli użyję sferycznego Mercatora, odległość kartezjańska wynosi 904 m, co jest daleka.

Tak, będziesz otrzymywać tego rodzaju błędy dzięki globalnej projekcji Mercatora: jest on dokładny na równiku, a zniekształcenie rośnie wykładniczo wraz z szerokością geograficzną od równika. Zniekształcenie odległości wynosi dokładnie 2 (100%) przy 60 stopniach szerokości geograficznej. Na waszych szerokościach testowych (64,14 stopni) obliczam zniekształcenie równe 2,294, dokładnie zgadzając się ze współczynnikiem 904/394 = 2,294. (Wcześniej obliczyłem 2,301, ale był on oparty na kuli, a nie na elipsoidzie WGS84. Różnica (0,3%) daje nam poczucie dokładności, jaką można uzyskać dzięki zastosowaniu projekcji opartej na elipsoidzie w porównaniu do opartej na kuli formuły Haversine. )

Nie ma czegoś takiego jak globalna projekcja, która zapewnia wszędzie bardzo dokładne odległości. To jeden z powodów, dla których używany jest system strefowy UTM!

Jednym z rozwiązań jest stosowanie geometrii sferycznej do wszystkich obliczeń, ale zostało to odrzucone (co jest rozsądne, jeśli zamierzasz wykonywać złożone operacje, ale warto ponownie podjąć decyzję).

Innym rozwiązaniem jest dostosowanie projekcji do porównywanych punktów . Na przykład możesz bezpiecznie użyć poprzecznego merkatora (jak w systemie UTM) z południkiem leżącym w pobliżu środka interesującego regionu. Przesuwanie południka jest proste: wystarczy odjąć długość południka od wszystkich długości i użyć pojedynczej projekcji TM wyśrodkowanej na południku głównym (ze współczynnikiem skali 1, a nie 0,9996 systemu UTM). W przypadku twojej pracy będzie to zwykle więcejdokładne niż przy użyciu samego UTM. Daje prawidłowe kąty (TM jest zgodny) i będzie niezwykle dokładne dla punktów oddzielonych zaledwie kilkadziesiąt kilometrów: spodziewaj się lepszej niż sześciocyfrowa dokładności. W rzeczywistości byłbym skłonny przypisywać wszelkie niewielkie różnice między tymi dostosowanymi odległościami TM a odległościami Haversine do różnicy między elipsoidą (używaną do projekcji TM) a kulą (używaną przez Haversine), zamiast do zniekształceń w występ.

Nie próbowałem tego, ale z dokumentacji wygląda na to, że można użyć http://search.cpan.org/~grahamc/Geo-Coordinates-UTM-0.08/UTM.pm#latlon_to_utm, aby uzyskać parę lat / lon (plus elipsoida) do strefy UTM i listy współrzędnych. Następnie możesz kontynuować obliczenia jak poprzednio.