Измерение расстояния в сферическом Меркаторе против зонного UTM

11

У меня есть очки в WGS84 широта / длина, и я хотел бы измерить «малые» (скажем, 5 км) расстояния между ними.

Я могу использовать формулу haversine из http://www.movable-type.co.uk/scripts/latlong.html, и она работает очень хорошо.

Я бы хотел использовать библиотеки Python Shapely, чтобы я мог выполнять больше операций, чем просто расстояние, и потому, что в масштабе, с которым я работаю, плоская земля является достаточно хорошим приближением. Чтобы надежно спроецировать географические координаты на декартову, я использую Python proj4, но, похоже, получаю больше ошибок, чем хотелось бы.

Если я использую локальную зону UTM, я получаю разницу между haversine в пару метров, и это нормально. Но я не хочу работать с зоной UTM (точки могут быть во всем мире), поэтому я попробовал использовать «сферический меркатор», но теперь разница между haversine и проецируемыми расстояниями составляет более 100%. Это действительно правильно для сферического Меркатора? Все, что я действительно хочу, это работоспособная декартова проекция для двух точек в пределах 5 км друг от друга в любой точке мира.

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)

В этот момент расстояние между ними составляет 394 м, а в зоне 27 utm - 395 м. Но если я использую сферический Меркатор, то декартово расстояние составляет 904 м, что далеко.

Карл П
источник
Зону UTM легко «проработать» на основе долгот. Выберите типичную долготу лямбда, -180 <= лямбда <180, и используйте ее для вычисления номера зоны как Int ((180 + лямбда) / 6) +1. Используйте знак широты, чтобы выбрать север и юг. Вам не нужно использовать специальные полярные зоны в высоких широтах; на самом деле, очень близко к полюсу вы можете использовать практически любую зону UTM.
whuber

Ответы:

17

Да, вы получите ошибки такого рода с глобальной проекцией Меркатора: она точна на экваторе, и искажение возрастает экспоненциально с широтой от экватора. Искажение расстояния составляет точно 2 (100%) при 60 градусах широты. На ваших тестовых широтах (64,14 градуса) я вычисляю искажение 2,294, точно согласуясь с соотношением 904/394 = 2,294. (Ранее я вычислял 2,301, но он основывался на сфере, а не на эллипсоиде WGS84. Разница (0,3%) дает нам ощущение точности, которую вы могли бы получить, используя проекцию на основе эллипсоидов по сравнению с формулой Хаверсайна на основе сфер. )

Не существует такой вещи, как глобальная проекция, которая повсюду дает очень точные расстояния. Это одна из причин использования системы зон UTM!

Одним из решений является использование сферической геометрии для всех ваших вычислений, но вы отклонили это (это разумно, если вы собираетесь выполнять сложные операции, но решение, возможно, стоит пересмотреть).

Другое решение - адаптировать проекцию к сравниваемым точкам . Например, вы могли бы безопасно использовать поперечный Меркатор (как в системе UTM) с меридианом, лежащим рядом с центром области интереса. Перемещение меридиана является простым делом: просто вычтите долготу меридиана из всех долгот и используйте одну проекцию TM, центрированную на главном меридиане (с масштабным коэффициентом 1, а не 0,9996 системы UTM). Для вашей работы это будет большеточнее, чем использование самого UTM. Это даст правильные углы (TM является конформным) и будет удивительно точным для точек, разделенных всего несколькими десятками километров: ожидайте лучше, чем шестизначная точность. Фактически, я был бы склонен приписывать любые небольшие различия между этими адаптированными расстояниями ТМ и расстояниями Хаверсайна различием между эллипсоидом (используемым для проекции ТМ) и сферой (используемым Хаверсайном), а не искажением в проекция.

Whuber
источник
Это звучит довольно прекрасно, я думаю, мне нужно создать собственную строку инициализации для proj4, а не использовать какие-либо из существующих строк EPSG?
Карл П
1
+1 адаптировать проекцию к точкам. Я предпочитаю поперечную пластину карре, а не поперечную проекцию Меркатора, но на достаточно небольших участках («крупномасштабный») почти любая проекция, «отцентрированная» возле интересующей области, даст хорошую точность.
Дэвид Кэри
@ Давид Интересная идея. На сфере поперечная пластина Карри (Кассини) будет близка к приближенной формуле, которую я дал на gis.stackexchange.com/posts/2964/edit (которая может быть приемлемым решением здесь). Формулы для TM и TPC похожи по сфере; на эллипсоиде TPC немного проще. ТМ, вероятно, поддерживается большим количеством программного обеспечения.
whuber
@Karl Вы можете использовать любую из зон TM, если хотите. Просто сдвиньте все долготы, чтобы центральная точка в интересующей вас области совпадала с центральным меридианом выбранной зоны. Умножьте все расстояния на 1 / 0,9996 (и умножьте все площади на квадрат этого фактора), не меняйте углы и ориентации и - если ваши вычисления дают новые точки - просто сдвиньте их долготы обратно к исходной системе координат ,
whuber
0

Я не пробовал этого, но из документации видно, что вы можете использовать http://search.cpan.org/~grahamc/Geo-Coordinates-UTM-0.08/UTM.pm#latlon_to_utm для получения пары широта / долгота (плюс эллипсоид) в зону UTM и список координат. Затем вы можете продолжить свой расчет, как и раньше.

Ян Тертон
источник