Если вы погружены в мир аналитики, скорее всего, вы сталкивались с разочаровывающими сценариями, когда расчет расстояния между двумя точками на Земле становится необходимостью. Однако вы быстро обнаруживаете серьезное препятствие — каждый ваш запрос имеет свою цену. Но не бойтесь, товарищи аналитики! В этом увлекательном блоге мы раскрываем подход, который может помочь вам рассчитать расстояние поездки с учетом маршрутов, ограничения скорости, одностороннего/двустороннего движения и других атрибутов, аналогичных картам Google, но бесплатно.
Случаи использования
Прежде чем погрузиться в технические подробности, давайте уделим немного времени тому, чтобы понять, в каких случаях полезно получить информацию о пройденном расстоянии.
- ETA (расчетное время прибытия) и автомобильные маршруты. В мире, где вы можете доставить почти все за считанные минуты, ETA (расчетное время прибытия) и автомобильные маршруты стали незаменимы для процветания бизнеса. Несомненно, использование платной услуги, такой как Google Distance Matrix API, кажется логичным выбором для бизнес-операций. Однако, когда дело доходит до проведения пилотных исследований или создания данных о расстоянии/маршруте исключительно для аналитических целей, инвестирование в дорогостоящие API может стать серьезной проблемой.
- Инженерия гиперлокальных особенностей.Представьте, что вы думаете об открытии ресторана в вашем районе, но не уверены в идеальном месте. В таком сценарии вам нужна исчерпывающая информация о потенциальном местоположении, чтобы принять взвешенное решение.
- Демографические данные
- Тип населения
- Уровень конкуренции
- Легкий доступ с автомагистралей, школ, жилых обществ и т. д.
Давайте сосредоточимся на одном факторе из этого примера: Конкуренция.
Традиционно вы можете нарисовать радиус в несколько километров от желаемого местоположения и подсчитать общее количество существующих ресторанов в этой области. Однако этот подход представляет собой проблему. Это может привести к переоценке конкуренции, потому что не все рестораны поблизости будут прямыми конкурентами.
Посмотрите на картинку ниже. Несмотря на то, что расстояние между двумя точками составляет всего 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 пунктов назначения, я также столкнулся с этой проблемой. Чтобы решить эту проблему, я разработал несколько приемов, которыми хочу поделиться с вами, чтобы сделать процесс быстрее и эффективнее.
Чтобы найти ближайшие узлы и рассчитать расстояния, нам понадобится фрейм геоданных или граф ребер. Есть два подхода, которые вы можете использовать:
- Загрузите весь фрейм геоданных (ребра) из osmnx для определенного города, над которым вы работаете. Кроме того, вы можете получить геоданные (ребра) для каждого пин-кода, сопоставив широту и долготу с соответствующим пин-кодом. Этот подход экономит время за счет создания единого фрейма геоданных (ребер) для нескольких строк, географически близких друг к другу.
- Используйте алгоритм кластеризации для исходных координат широты и долготы. Это позволяет группировать источники, которые находятся в непосредственной близости друг от друга, в кластеры. Затем вы можете создать один кадр геоданных (ребра) для всех источников в каждом кластере.
- Не стройте каждый из графиков при выполнении расчетов.
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. Я более чем счастлив помочь вам и предоставить дальнейшие рекомендации. Мы очень ценим ваши отзывы и участие, поскольку они способствуют созданию полезного и поддерживающего сообщества.