Вычисление средней / длинной точки Python дает неверный результат, когда долгота> 90

Question

У меня проблема с короткой функцией для вычисления средней точки линии при заданной широте и долготе точек на каждом конце. Проще говоря, он работает правильно, когда долгота больше -90 градусов или меньше 90 градусов. Для другой половины планеты это дает несколько случайный результат.

Код представляет собой Python-преобразование javascript, предоставленного в http://www.movable -type.co.uk / scripts / latlong.html и, похоже, соответствует исправленным версиям здесь и здесь . Сравнивая с двумя версиями stackoverflow, я признаю, что не кодирую на C # или Java, но не могу определить, где моя ошибка.

Код выглядит следующим образом:

#!/usr/bin/python

import math

def midpoint(p1, p2):
   lat1, lat2 = math.radians(p1[0]), math.radians(p2[0])
   lon1, lon2 = math.radians(p1[1]), math.radians(p2[1])
   dlon = lon2 - lon1
   dx = math.cos(lat2) * math.cos(dlon)
   dy = math.cos(lat2) * math.sin(dlon)
   lat3 = math.atan2(math.sin(lat1) + math.sin(lat2), math.sqrt((math.cos(lat1) + dx) * (math.cos(lat1) + dx) + dy * dy))
   lon3 = lon1 + math.atan2(dy, math.cos(lat1) + dx)
   return(math.degrees(lat3), math.degrees(lon3))

p1 = (6.4, 45)
p2 = (7.3, 43.5)
print "Correct:", midpoint(p1, p2)

p1 = (95.5,41.4)
p2 = (96.3,41.8)
print "Wrong:", midpoint(p1, p2)

Есть предложения?

Answer 1

Замените ваш код настройки arg на:

lat1, lon1 = p1
lat2, lon2 = p2
assert -90 <= lat1 <= 90
assert -90 <= lat2 <= 90
assert -180 <= lon1 <= 180
assert -180 <= lon2 <= 180
lat1, lon1, lat2, lon2 = map(math.radians, (lat1, lon1, lat2, lon2))

и снова запустите ваш код.

Обновите Несколько общих рекомендаций о вычислениях, связанных с широтой/ долгота:

1. Входной широта / долгота в градусах или радианах?
2. Проверить входной широта / долгота для допустимого диапазона
3. Проверить ВЫХОДНОЙ широта / долгота для действительного диапазона.Долгота имеет разрыв на международной линии времени.

Последнюю часть процедуры средней точки можно было бы с пользой изменить, чтобы избежать потенциальной проблемы при использовании на дальние расстояния:

lon3 = lon1 + math.atan2(dy, math.cos(lat1) + dx)
# replacement code follows:
lon3d = math.degrees(lon3)
if lon3d < -180:
    print "oops1", lon3d
    lon3d += 360
elif lon3d > 180:
    print "oops2", lon3d
    lon3d -= 360
return(math.degrees(lat3), lon3d)

ДляНапример, нахождение средней точки между Оклендом, Новая Зеландия (-36,9, 174,8) и Папеэте, Таити (-17,5, -149,5) дает oops2 194.270430902 на пути к правильному ответу (-28.355951246746923, -165.72956909809082)

Answer 2

Прежде всего, извиняюсь, что оставляю другой ответ.У меня есть решение проблемы, упомянутой в этом ответе, заключающееся в том, как найти среднюю точку двух точек по обе стороны от линии даты.Я хотел бы просто добавить комментарий к существующему ответу, но у меня нет репутации, чтобы сделать это.

Решение было найдено, посмотрев на файлы Javascript, включающие инструмент на http://www.movable -type.co.uk / скрипты / latlong.html .Я использую shapely.geometry.Point.Если вы не хотите устанавливать этот пакет, то использование кортежей вместо этого будет работать точно так же.

def midpoint(pointA, pointB):
    lonA = math.radians(pointA.x)
    lonB = math.radians(pointB.x)
    latA = math.radians(pointA.y)
    latB = math.radians(pointB.y)

    dLon = lonB - lonA

    Bx = math.cos(latB) * math.cos(dLon)
    By = math.cos(latB) * math.sin(dLon)

    latC = math.atan2(math.sin(latA) + math.sin(latB),
                  math.sqrt((math.cos(latA) + Bx) * (math.cos(latA) + Bx) + By * By))
    lonC = lonA + math.atan2(By, math.cos(latA) + Bx)
    lonC = (lonC + 3 * math.pi) % (2 * math.pi) - math.pi

    return Point(math.degrees(lonC), math.degrees(latC))

Я надеюсь, что это полезно и не считается неуместным, поскольку это ответ на вопрос.поднят в предыдущем ответе.

Answer 3

Не прямой ответ на вопрос> 90, но это помогло в моем случае, поэтому я оставляю это здесь на всякий случай, если это может помочь другим.

Мне нужно было только поставить серединууказать на карте лат, длинный в градусах.Я не использовал преобразование в радианы, и оно правильно отображено на карте с использованием простого евклидова расстояния.Пример функции:

def midpoint_euclidean(x1,y1,x2,y2):
    dist_x = abs(x1-x2) / 2.
    dist_y = abs(y1-y2) / 2.
    res_x = x1 - dist_x if x1 > x2 else x2 - dist_x
    res_y = y1 - dist_y if y1 > y2 else y2 - dist_y
    return res_x, res_y