Если вы погружены в мир аналитики, скорее всего, вы сталкивались с разочаровывающими сценариями, когда расчет расстояния между двумя точками на Земле становится необходимостью. Однако вы быстро обнаруживаете серьезное препятствие — каждый ваш запрос имеет свою цену. Но не бойтесь, товарищи аналитики! В этом увлекательном блоге мы раскрываем подход, который может помочь вам рассчитать расстояние поездки с учетом маршрутов, ограничения скорости, одностороннего/двустороннего движения и других атрибутов, аналогичных картам Google, но бесплатно.

Случаи использования

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

  1. ETA (расчетное время прибытия) и автомобильные маршруты. В мире, где вы можете доставить почти все за считанные минуты, ETA (расчетное время прибытия) и автомобильные маршруты стали незаменимы для процветания бизнеса. Несомненно, использование платной услуги, такой как Google Distance Matrix API, кажется логичным выбором для бизнес-операций. Однако, когда дело доходит до проведения пилотных исследований или создания данных о расстоянии/маршруте исключительно для аналитических целей, инвестирование в дорогостоящие API может стать серьезной проблемой.
  2. Инженерия гиперлокальных особенностей.Представьте, что вы думаете об открытии ресторана в вашем районе, но не уверены в идеальном месте. В таком сценарии вам нужна исчерпывающая информация о потенциальном местоположении, чтобы принять взвешенное решение.
  • Демографические данные
  • Тип населения
  • Уровень конкуренции
  • Легкий доступ с автомагистралей, школ, жилых обществ и т. д.

Давайте сосредоточимся на одном факторе из этого примера: Конкуренция.
Традиционно вы можете нарисовать радиус в несколько километров от желаемого местоположения и подсчитать общее количество существующих ресторанов в этой области. Однако этот подход представляет собой проблему. Это может привести к переоценке конкуренции, потому что не все рестораны поблизости будут прямыми конкурентами.
Посмотрите на картинку ниже. Несмотря на то, что расстояние между двумя точками составляет всего 200 метров, пройденное расстояние составляет 2200 метров. В подобных случаях для точной оценки конкуренции более целесообразно учитывать дальность поездки.
Не только для конкуренции, но и для четвертой точки, которая является легкодоступной. пройденное расстояние» имеет больше смысла, чем кратчайшее расстояние.

Пройденное расстояние

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

import numpy as np
import pandas as pd
import geopandas as gpd
from shapely.geometry import LineString, Point
import osmnx  as ox
import networkx as nx
import matplotlib.pyplot as plt
import pickle
import warnings
warnings.filterwarnings("ignore")

Начнем с выбора двух точек на карте:

Начальная точка: Domino's Pizza в торговом центре Vega City Mall, Бангалор, Индия (12,9073 с.ш., 77,6011 в.д.)
Пункт назначения: Domino's Pizza на 2-м этапе BTM, Бангалор, Индия (12,9197 с.ш., 77,6125 в.д.)

Прежде чем углубиться в расстояние пути, мы сначала рассчитаем расстояние по гаверсинусу. Этот метод широко используется из-за его простоты и легкости расчета в случаях, подобных нашему.

def haversine_np(lon1, lat1, lon2, lat2):
    
    """
    Calculate the great circle distance between two points
    on the earth (specified in decimal degrees)

    All args must be of equal length.    

    """
    lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])

    dlon = lon2 - lon1
    dlat = lat2 - lat1

    a = np.sin(dlat/2.0)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2.0)**2

    c = 2 * np.arcsin(np.sqrt(a))
    d = (6371 * c)
    return round (d*1000,3)

g1=(77.6011,12.9073)
g2=(77.6125,12.9197)

haversine_np(g1[1],g1[0],g2[1],g2[0])

# Distance is 1301.705 meters

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

Теперь давайте рассчитаем расстояние перемещения и сравним его с расстоянием гаверсинуса. Для начала нам нужно создать кадры геоданных для исходной и конечной точек. Кадр геоданных — это, по сути, кадр данных pandas с дополнительным столбцом под названием «геометрия». Столбец геометрии представляет форму каждой строки, которая может быть точкой (определенное место на Земле), линии (например, дороги или пляжи), или многоугольник (представляющий города, штаты, страны и т. д.). В нашем случае мы должны преобразовать координаты источника и пункта назначения в геометрию точки, чтобы рассчитать расстояние перемещения.

def lat_long_to_gdf(g1,g2):

    #Origin
    origin = gpd.GeoDataFrame(columns = ['name', 'geometry'], crs = 4326, geometry = 'geometry')
    origin.at[0, 'geometry'] = Point(g1[1] , g1[0])
    origin.at[0, 'name'] = 'Origin'
    # dest
    destination = gpd.GeoDataFrame(columns = ['name', 'geometry'], crs = 4326, geometry = 'geometry')
    destination.at[0, 'geometry'] = Point(g2[1],g2[0])
    destination.at[0, 'name'] = 'Destination'
    
    return origin,destination

origin,destination = lat_long_to_gdf(g1,g2)
print(origin)
print(destination)

Получив два кадра геоданных, давайте соединим их оба, чтобы создать полигон, соединяющий две точки с буфером 1000 метров.

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

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

График

Фрейм геоданных «Ребра», полученный из osmnx, содержит обширную информацию об области. Как упоминалось ранее, он состоит из нескольких узлов, распределенных по картируемой области, причем каждый узел связан с определенными координатами широты и долготы (представленными столбцами «u» и «v» в кадре геоданных ребер).

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

Гаверсинус против расстояния перемещения

Верхнее расстояние между двумя точками составляло 1300 метров, но расстояние в пути составляет 2300 метров приблизительно.
Это ясно показывает, почему и как верное расстояние может вводить в заблуждение. много раз

Функция поиска кратчайшего пути и построения карты

def shortest_path_viz2(origin, destination,buf = 0.01 network = 'drive'):
    '''
    origin and destination  
    crs 4326, 
    buf 0.01 (1000 meters buffer around the origin and des points)
    netwrok string {"all_private", "all", "bike", "drive", "drive_service", "walk"}
    
    return shortest path 
    '''
    len_list = []
    
    # creating a frame for fetching OSM data
    combined = origin.append(destination)    
    convex = combined.unary_union.convex_hull        
    graph_extent = convex.buffer(buf)
    graph = ox.graph_from_polygon(graph_extent, network_type= network)
    
    # Reproject the graph
    graph_proj = ox.project_graph(graph)    
    # Get the GeoDataFrame
    edges = ox.graph_to_gdfs(graph_proj, nodes=False)   
    # Get CRS info UTM
    CRS = edges.crs
    
    # Reproject all data
    origin_proj = origin.to_crs(crs=CRS)
    destination_proj = destination.to_crs(crs=CRS)
    
    # routes of shortest path
    routes = gpd.GeoDataFrame()
    
    # Get nodes from the graph
    nodes = ox.graph_to_gdfs(graph_proj, edges=False)
    print("third")
    
    # Iterate over origins and destinations
    for oidx, orig in origin_proj.iterrows():
        
        # Find closest node from the graph --> point = (latitude, longitude)
        closest_origin_node = ox.distance.nearest_nodes(G=graph_proj, X = orig.geometry.x, Y = orig.geometry.y)
        
        # Iterate over targets
        for tidx, target in destination_proj.iterrows():
            # Find closest node from the graph --> point = (latitude, longitude)
            closest_target_node = ox.distance.nearest_nodes(graph_proj, Y = target.geometry.y, X= target.geometry.x)
            
            try:
                # Check if origin and target nodes are the same --> if they are --> skip
                if closest_origin_node == closest_target_node:
                    print("Same origin and destination node. Skipping ..")
                    continue

                # Find the shortest path between the points
                route = nx.shortest_path(graph_proj, 
                                         source=closest_origin_node, 
                                         target=closest_target_node, weight='length')
                #length
                length = nx.shortest_path_length(graph_proj, source=closest_origin_node, target=closest_target_node, weight='length')
                len_list.append(length)
                # Extract the nodes of the route
                route_nodes = nodes.loc[route]

                # Create a LineString out of the route
                path = LineString(list(route_nodes.geometry.values))

                # Append the result into the GeoDataFrame
                routes = routes.append([[path]], ignore_index=True)
            except:
                continue

    # Add a column name
    routes.columns = ['geometry']        
    # Set geometry
    routes = routes.set_geometry('geometry')    
    # Set coordinate reference system
    routes.crs = nodes.crs    
    # Add a column name
    routes.columns = ['geometry']
    # Set geometry
    routes = routes.set_geometry('geometry')
    # Set coordinate reference system
    routes.crs = nodes.crs   
    plt.style.use('seaborn')
    # Plot
    ax = edges.plot(figsize=(16, 10), color='gray', linewidth=0.5, alpha=0.7)
    ax = origin_proj.plot(ax=ax, color='red')
    ax = destination_proj.plot(ax=ax, color='blue')
    ax = routes.plot(ax=ax, linewidth=3, alpha = 0.8, color = 'magenta')
    plt.axis('off')
        
    return CRS,graph_proj

Важные вещи, которые нужно помнить

Если вы до сих пор следили за этим блогом, вы могли заметить недостаток обсуждаемого алгоритма, а именно его скорость.

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

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

  1. Загрузите весь фрейм геоданных (ребра) из osmnx для определенного города, над которым вы работаете. Кроме того, вы можете получить геоданные (ребра) для каждого пин-кода, сопоставив широту и долготу с соответствующим пин-кодом. Этот подход экономит время за счет создания единого фрейма геоданных (ребер) для нескольких строк, географически близких друг к другу.
  2. Используйте алгоритм кластеризации для исходных координат широты и долготы. Это позволяет группировать источники, которые находятся в непосредственной близости друг от друга, в кластеры. Затем вы можете создать один кадр геоданных (ребра) для всех источников в каждом кластере.
  3. Не стройте каждый из графиков при выполнении расчетов.
graph = ox.graph_from_place('Bangalore', network_type='drive')

Давайте рассмотрим пример, в котором у вас есть 1000 источников, и вы можете объединить их в 10 кластеров. Если каждый кластер покрывает площадь 10 000 метров, вы можете воспользоваться преимуществами этой кластеризации для оптимизации своих вычислений.

Вместо создания отдельных краевых фреймов геоданных для всех 1000 источников вам нужно создать всего 10 краевых фреймов геоданных, каждый из которых представляет кластер источников в указанной области. Делая это, вы значительно снижаете вычислительную нагрузку и экономите время, так как можете выполнять вычисления для каждого кластера в целом, а не для отдельных источников. Такой подход помогает повысить эффективность и скорость обработки при работе с большим количеством источников или пунктов назначения.

Оптимизированный алгоритм

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

import time
import itertools
from collections import Counter

start_time = time.time()

def shortest_path_opt2(origin, destination,road = False, network = 'drive'):

        skip = False
        node_dist = {}
        
        try:
            
            start_time = time.time()
            # creating a frame for fetching OSM data
            combined = origin.append(destination)

            convex = combined.unary_union.convex_hull

            graph_extent = convex.buffer(0.03)

            # fetching graph
            graph = ox.graph_from_polygon(graph_extent, network_type= network)

            # Reproject the graph
            graph_proj = ox.project_graph(graph)
            
            # Get the GeoDataFrame
            edges = ox.graph_to_gdfs(graph_proj, nodes=False)

            # Get CRS info UTM
            CRS = edges.crs

            # Reproject all data
            origin_proj = origin.to_crs(crs=CRS)
            destination_proj = destination.to_crs(crs=CRS)
            
            end_time = time.time()
            elapsed_time = end_time - start_time
            print(f"time taken in df preprocess for {len(destination)} rows is : ", elapsed_time)
            
            start_time = time.time()
            try:
                
                # Find closest node from the graph --> point = (latitude, longitude)
                closest_origin_node = ox.distance.nearest_nodes(graph_proj, X = origin_proj.geometry.x, Y = origin_proj.geometry.y)
                
            except Exception as e:
                print('no origin node')
                print(e)
                destination['travel_dist'] = ["no origin node"]*len(destination)
                return destination
                
            try:
                
                # Find closest node from the graph --> point = (latitude, longitude)
                closest_target_node = ox.distance.nearest_nodes(graph_proj, X= destination_proj.geometry.x, Y = destination_proj.geometry.y)
                end_time = time.time()
                elapsed_time = end_time - start_time
                print(f"time taken in finding nearest node is : ", elapsed_time)

                destination['nodes'] = closest_target_node
                closest_target_node_unique = set(closest_target_node)


            except Exception as e:
                print('no dest node')
                print(e)
                destination['travel_dist'] = ["no dest node"]*len(destination)
                return destination

            start_time = time.time()
            
            #shortest_path_lengths = nx.multi_target_dijkstra_path_length(graph_proj, origin_nodes, target_node, weight='length')
            #shortest_path_lengths = nx.multi_source_dijkstra_path_length(graph_proj, closest_target_node, target_node, weight='length')

            for i in closest_target_node_unique:

                try:
                    
                    # Check if origin and target nodes are the same --> if they are --> skip
                    if closest_origin_node[0] == i:
                        print("Same origin and destination node. Skipping ..")
                        node_dist[i] = 0
                        continue
                        
                    if road:
                        length = nx.shortest_path_length(graph_proj, source=i,target=closest_origin_node[0], weight='length')
                    else:
                        length = nx.shortest_path_length(graph_proj, source=closest_origin_node[0],target=i, weight='length')

                    node_dist[i] = length
                    
                except Exception as e:
                    print('no route')
                    print(e)
                    node_dist[i] = 'no route'

            end_time = time.time()
            elapsed_time = end_time - start_time
            print(f"time taken in finding shortest route is : ", elapsed_time)
            
            
            start_time = time.time()
            
            temp = pd.DataFrame(node_dist.items(), columns=['nodes', 'travel_dist'])            
            destination = destination.merge(temp, on= 'nodes', how = 'left')
            
            end_time = time.time()
            elapsed_time = end_time - start_time
            print(f"time taken in final merge : ", elapsed_time)
            
            
        except Exception as e:
            print(e)
            destination['travel_dist'] = ["Graph not found"]*len(destination)


        return destination

Заключение

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

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