«Могу делать сайты, а могу и не делать»
Запись от 2 июня 2017 г.
Перевод географических координат в расстояние и обратно на Python

Итак, задача: есть куча географических координат объектов, необходимо перевести эти координаты в расстояние от какой-то определенной точки. Зачем это может понадобиться? Например, вы хотите разбить какие-то объекты в области на зоны, проведя кластеризацию облака точек.

Что мы знаем о земле? Земля есть шар, значит она обладает определенным радиусом, который в среднем равен 6371 километр. Кроме радиуса земли понадобятся координаты начальной точки. Это может быть как определенная точка в облаке, так и произвольная.

Алгоритм вычисление расстояния от начальной точки до искомой достаточно прост: вычисляем сначала расстояние Х, изменяя только долготу, затем вычисляем расстояние Y, изменяя только широту. Для вычисления расстояния воспользуемся методом вычисления длины дуги большого круга. Если верить географии, между двумя точками, если они не авляются антиподами друг друга, можно провести уникальный большой круг, а потом измерить длину кратчайшей дуги между искомыми точками. Это сделать можно разными способами. Например, с помощью гаверсинусов.

Реализация формулы гаверсинусов на python:

def haversine(lon1, lat1, lon2, lat2):
    """
    Calculate the great circle distance between two points
    on the earth (specified in decimal degrees)
    :param lon1: first point longitude
    :param lat1: first point latitude
    :param lon2: second point longitude
    :param lat2: second point latitude
    :return: distance between two points in meters
    """
    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
    c = 2 * asin(sqrt(a))
    dec_m = EARTH_RADIUS * c
    return dec_m

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

Расстояние находить научились, теперь можно перейти к вычислениям координат точки:

def translate_geo_to_decart(data):
    """
    converts lonlat into meters
    :param data: dict of data in format [[lon, lat], ...]
    :return: converted data in format [[x, y], ...]
    """
    buf = []
    for point in data:
        x = haversine(float(START_POINT[1]), float(START_POINT[0]),
                      float(START_POINT[1]), float(point[0]))
        y = haversine(float(START_POINT[1]), float(START_POINT[0]),
                      float(point[1]), float(START_POINT[0]))
        buf.append([x, y])
    return buf

Функция получает на вход массив точек с географическими координатами и вычисляет расстояние по X и Y каждой от начальной точки. Но раз прямой перевод делать умеем, то необходимо уметь делать и обратный, ведь так? Алгоритм работы анологичный: вычисляем длину дуги в метрах, проводим большой круг и обратным гаверсинусом вычисляем разницу широт и долгот для искомой точки. Прибавляя полученные значения к координатам начальной точки, получаем искомую координату:

def translate_geo_to_decart(data):
    """
    converts lonlat into meters
    :param data: dict of data in format [[lon, lat], ...]
    :return: converted data in format [[x, y], ...]
    """
    buf = []
    for point in data:
        x = haversine(float(START_POINT[1]), float(START_POINT[0]),
                      float(START_POINT[1]), float(point[0]))
        y = haversine(float(START_POINT[1]), float(START_POINT[0]),
                      float(point[1]), float(START_POINT[0]))
        buf.append([x, y])
    return buf

Полный текст программы:

from math import radians, sin, cos, asin, sqrt, degrees


START_POINT = [54.24, 35.13]
EARTH_RADIUS = 6371 * 100


def haversine(lon1, lat1, lon2, lat2):
    """
    Calculate the great circle distance between two points
    on the earth (specified in decimal degrees)
    :param lon1: first point longitude
    :param lat1: first point latitude
    :param lon2: second point longitude
    :param lat2: second point latitude
    :return: distance between two points in meters
    """
    lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
    c = 2 * asin(sqrt(a))
    dec_m = EARTH_RADIUS * c
    return dec_m


def translate_geo_to_decart(data):
    """
    converts lonlat into meters
    :param data: dict of data in format [[lon, lat], ...]
    :return: converted data in format [[x, y], ...]
    """
    buf = []
    for point in data:
        x = haversine(float(START_POINT[1]), float(START_POINT[0]),
                      float(START_POINT[1]), float(point[0]))
        y = haversine(float(START_POINT[1]), float(START_POINT[0]),
                      float(point[1]), float(START_POINT[0]))
        buf.append([x, y])
    return buf


def reverse_transofm_coordinates(data):
    """
    reverse translate coords into [lat, lon] from distances
    :param data: dict of data in format [[x, y], ...]
    :return: point in format [lat lon]
    """
    buf = []
    for point in data:
        distance = sqrt(float(point[0])**2 + float(point[1])**2)
        dlat = distance / EARTH_RADIUS
        dlon = asin(sin(dlat) / cos(radians(START_POINT[1])))
        buf.append([START_POINT[0] + degrees(dlon), START_POINT[1] + degrees(dlat)])
    return buf

Ссылка на Github: https://github.com/whspr/geographic-coords-transform