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

11

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.

Karl P.
źródło
Strefę UTM można łatwo „opracować” w oparciu o długości geograficzne. Wybierz typową lambda długości geograficznej, -180 <= lambda <180 i użyj jej do obliczenia numeru strefy jako Int ((180 + lambda) / 6) +1. Użyj znaku szerokości geograficznej, aby wybrać między północą a południem. Nie musisz używać specjalnych stref polarnych na dużych szerokościach geograficznych; tak naprawdę bardzo blisko bieguna można użyć prawie każdej strefy UTM.
whuber

Odpowiedzi:

17

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.

Whuber
źródło
Brzmi to całkiem perfekcyjnie, chyba muszę stworzyć własny ciąg inicjujący dla proj4, zamiast móc korzystać z istniejących ciągów EPSG?
Karl P
1
+1 dostosuj rzut do punktów. Wolę poprzeczną płytę Carrée niż poprzeczną Mercator, ale w przypadku wystarczająco małych obszarów („duża skala”), prawie każdy występ „wyśrodkowany” w pobliżu obszaru zainteresowania zapewni dobrą dokładność.
David Cary,
@David Ciekawy pomysł. Na kuli poprzeczna płytka Carree (Cassini) będzie zbliżona do przybliżonej formuły, którą podałem na gis.stackexchange.com/posts/2964/edit (co może być tutaj akceptowalnym rozwiązaniem). Formuły dla TM i TPC są podobne na kuli; na elipsoidzie TPC jest nieco prostsze. TM jest prawdopodobnie obsługiwany przez więcej programów.
whuber
@Karl Możesz użyć dowolnej strefy TM, jeśli chcesz. Po prostu przesuń wszystkie długości geograficzne, aby centralny punkt w twoim regionie zainteresowania pokrywał się z centralnym południkiem wybranej strefy. Pomnóż wszystkie odległości przez 1 / 0,9996 (i pomnóż wszystkie obszary przez kwadrat tego współczynnika), nie zmieniaj żadnych kątów ani łożysk, i - jeśli twoje obliczenia dają nowe punkty - po prostu przesuń ich długości z powrotem do pierwotnego układu współrzędnych .
whuber