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

Часть первая. Два способа найти на карте расстояние от точки до области

Если ваше мобильное приложение действительно мобильное, оно работает с координатами устройства. Местоположение пользователя (и устройства) влияет на различные бизнес-показатели приложения, такие, как стоимость доставки, коэффициент сложности работ и т.п.
Ниже я покажу примеры реализации алгоритмов на python с использованием библиотек scipy и shapely.
Для геокодирования мы используем Google Maps. Он подходит и по функционалу, и по стоимости использования. На момент написания этой статьи Google разрешает бесплатно сделать первые 20 000 запросов геокодинга в месяц.

Способ 1. Рассчитываем маршрут на основании вершин полигона

Предположим, нам надо найти расстояние от точки где-то в Московской области до МКАД. Нужен реальный путь, а не геометрический. Поэтому сначала строим полигон из точек выездов с МКАД, а они не совпадают с вершинами очертания дороги на карте.

exit_coordinates: List[Tuple[float, float]]
latitude: float
longitude: float

Для работы с геометрией используем библиотеку shapely. С её помощью легко определить, находится ли интересующая нас точка внутри полигона, или нет. Если находится — расстояние очевидно равно 0.

from shapely.geometry import Polygon, Point
polygon = Polygon(exit_coordinates)
polygon.contains(Point(latitude,longitude))

Нас интересует случай, когда точка находится вне полигона. Задача найти ближайшие объекты (выезды с МКАД) до исходной точки — довольно типичная в математике. Обычно её решают с помощью KD-дерева. Так и поступим. В питоне дерево реализовано в библиотеке scipy. Найдём ближайшие вершины из выездов с МКАД до нашей точки.

from scipy.spatial import KDTree
kd_tree = KDTree(exits_coordinates)
dists, indexes = kd_tree.query((latitude, longitude), k=3)
nearest_coordinates = list()
for _, index in zip(dists, indexes):
		nearest_coordinates.append(exits_coordinates[index])

Число k — количество вершин — не должно быть слишком маленьким, ведь ближайший выезд со МКАД может быть временно перекрыт. Или может находиться где-нибудь за рекой и не иметь прямого пути до нашей точки. Слишком большое k тоже для нас бесполезно, т.к. подходящий выезд находится в нескольких ближайших вершинах. Теперь обратимся к сервису Google Maps. В нём уже есть функционал поиска расстояний из массива точек до конкретной точки — Distance Matrix API.
Upd: сервис Distance Matrix API биллингует маршрут до каждой ближайшей вершины отдельно. Таким образом, этот способ обходится дороже, чем второй, описанный далее. Спасибо sovetnikov

import googlemaps
gmaps_client = googlemaps.Client(key=settings.GOOGLE_MAPS_API_KEY)
gmaps_client.distance_matrix(
    origins=nearest_coordinates,
    destinations=(latitude, longitude),
    mode='driving',
)

Нам остаётся только распарсить ответ. Обычно в нём несколько маршрутов, и самый первый не всегда самый короткий. Поэтому сами выбираем подходящее значение.

Способ 2. Рассчитываем маршрут от центра полигона

Теперь помимо вершин полигона мы определяем некоторую точку внутри него, из которой будем строить все маршруты Google Maps. Воспользуемся Directions API, который построит для нас маршруты, и выберем кратчайший из них.

start_point: Tuple[float, float],
destination: Tuple[float, float]
polygon: shapely.geometry.Polygon
gmaps_client = googlemaps.Client(key=settings.GOOGLE_MAPS_API_KEY)
driving_path = gmaps_client.directions(start_point, destination)

Начинаем итеративно с конца складывать длины отрезков пути то тех пор, пока координата начала отрезка не окажется внутри полигона (парсинг опущен):

for step in reversed(driving_path):
    start_point = step['start_location']['lat'], step['start_location']['lng']
    if is_inside_polygon(start_point, self.polygon):
        end_point = step['end_location']['lat'], step['end_location']['lng']
        distance += get_part_outside_polygon(
            get_line(start_point, end_point), polygon
        ) * step['distance']['value']
        break
    distance += step['distance']['value']

Нам остается только прибавить часть отрезка вне полигона. Для этого с помощью библиотеки shapely строим геометрическую линию между координатами начала и конца пути и находим, какая часть её длины лежит вне полигона. Такой же процент от длины высчитываем для полученного отрезка пути.

Путь — это маршрут на местности по реальным дорогам с естественной кривизной. Если это длинная прямая (проспект или магистраль), то наше приближение позволит высчитать маршрут точно в процентном соотношении.
Если же пересекающая полигон дорога достаточно кривая, такое приближение будет неверным. Но искривлённые части в самом маршруте Google Maps обычно имеют малую длину, и погрешности их вычисления не повлияют результат.

from shapely.geometry import LineString, Polygon
line = LineString(point1, point2)
part_outside_len = float(line.difference(polygon).length) / float(line.length)

Если честно, я не сравнивал эти два метода. Пользуюсь ими больше года, оба работают без сбоев. Решил не расписывать подробно их реализацию. Вместо этого открыл доступ к своей библиотеке по работе с гео. Либа умеет и в обычный геокодинг, в том числе с эффективным кешированием. Issues и pull-requests приветствуются!

Часть вторая. Экономим запросы обратного геокодирования

Часто нам нужно определить адреса точек, находящиеся рядом и соответствующие одному объекту в реальном мире. Каждый раз делать запрос во внешнюю систему геокодирования не круто, и тут резонно возникает вопрос о кешировании.
Обычно клиент присылает координаты с избыточной точностью, например, 59.80447691, 30.39570337. Но сколько знаков в дробной части будут для нас значимыми?
Для ленивых и торопящихся ответ — четыре. Для всех остальных — пояснение ниже.

Сначала немного географии.

• Длина экватора — 40 075,696 км, и это нулевая широта.
• Меридианы — линии долготы, перпендикулярны линиям широты. Длина любого меридиана — 40 008,55 км
• Градус широты — 40 008,55 км / 360 = 111,134861111 км. А значит, одна сотая дает изменение приблизительно в километр, а одна десятитысячная — в 10 метров.
• Длина окружности линии широты уменьшается от экватора, и находится умножением на косинус угла широты, следовательно, для 60 градусов широты (широты Санкт-Петербурга) длина меньше ровно в 2 раза. Тогда градус долготы — 40 075,696 * 0.5 / 360 = 55,66066888 км, или одна десятитысячная — 5 метров.

Погрешность в 1/10000 градуса дает прямоугольник погрешности 5-10 на 10 метров. Это позволит нам точно «попасть» в здание, т.к. здание значительно крупнее точки. А если из-за погрешности точка не попадает в здание, Google всё равно определяет ближайшее к ней.
Округлять до трёх знаков рискованно, точность снижается уже существенно. А до пяти — не имеет смысла, потому что точность GPS в смартфонах обычно не превышает всего нескольких метров.

Таким образом, если в кеше есть значение с ключом «address:59.8045,30.3957», то все координаты, при округлении которых до 4х знаков получается такой же ключ, соответствуют одному геокодированному адресу. Делаем меньше запросов — работаем быстрее и меньше платим за использование сервиса геокодинга. Профит!