Планируем идеальный поход с NetworkX и OpenStreetMap

Любой заядлый турист знает: всё, что ждёт его в походе, зависит от подготовки. Брать ли дождевики? Сколько идти от одного источника воды до другого? Где лучше ночевать в этой местности в это время года? Но самый важный вопрос звучит проще простого: «Что я там буду делать?» Отвечаем на эти вопросы к старту нашего курса по Fullstack-разработке на Python.

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

Поэтому я решил создать собственное приложение для планирования походов с открытым исходным кодом и механизмом оптимизации маршрута, которое дополнит функционал коммерческих сайтов. В этой статье я проложу (простите за каламбур) кратчайший маршрут по всем независимым компонентам системы и покажу их в совокупности в демоверсии. Хотите увидеть код? Добро пожаловать в реп — вот сюда.

Когда нужны бесплатные данные для выбора маршрута, очевидный выбор — карты OpenStreetMap. Даже AllTrails берёт именно их базу дорог через Overpass API и отсеивает автомобильные дороги, оставляя туристические тропы.

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

Как и все дороги, туристические тропы имеют перекрёстки (вершины графа) и отрезки пути (рёбра графа). А значит, для выбора оптимального маршрута разумно использовать метод графов, например алгоритм Дейкстры.

Америку я вам тут не открою, ведь обучающих материалов и библиотек по построению маршрута на графах создано великое множество. Такие пакеты, как OSMnx, позволяют формировать высокоуровневые запросы к наборам данных OSM по геокоду (например, «Shenandoah National Park, Virginia, USA») и конвертировать результат в графы NetworkX.

Вот листинг функции load_map, с помощью которой я получил данные о тропах и сохранял их в виде графов (для эффективного расчёта) и файлов GeoJSON (для визуализации).

def load_map(place: str, graph_only: bool = False) -> (nx.MultiDiGraph, dict):     """     Load OSM trail data for a given region. Initially check if the graph has already been cached on disk, otherwise it will be downloaded.     Args:         place: The geocode of the region of interest, e.g. 'Shenandoah National Park, Virginia, USA'         graph_only: If true, return only the NetworkX graph, not the geojson.      Returns: The dataset as a NetworkX MultiGraph, the nodes geojson, the edges geojson      """     filename = p.sub("", place).lower().replace(" ", "_")     graph_fp = f"static/graphs/{filename}.graphml"     edges_fp = f"static/edges/{filename}.geojson"     nodes_fp = f"static/nodes/{filename}.geojson"     try:         G = ox.load_graphml(graph_fp)         print("Graph loaded from disk")     except FileNotFoundError:         print("Downloading and caching graph from OSM")         # custom filter to only include walkable segments - see: https://support.alltrails.com/hc/en-us/articles/360019246411-OSM-Derivative-Database-Derivation-Methodology         cf = '["highway"~"path|track|footway|steps|bridleway|cycleway"]'         G = ox.graph_from_place(place, custom_filter=cf)         ox.save_graphml(G, graph_fp)     if graph_only:         return G     if os.path.isfile(edges_fp) is False or os.path.isfile(nodes_fp) is False:         # convert Graph to GeoDataFrames and save as GeoJSON         nodes, edges = ox.graph_to_gdfs(G)         edges = edges[["name", "geometry", "length"]]         edges["name"] = edges["name"].apply(lambda x: x[0] if type(x) is list else x)         edges.to_file(edges_fp)         nodes.to_file(nodes_fp)     with open(nodes_fp, "r") as f:         nodes = json.load(f)     with open(edges_fp, "r") as f:         edges = json.load(f)     for feat in edges["features"]:         # add custom tooltip property for Leaflet visualization         tooltip = f"{feat['properties']['name']}, {round(feat['properties']['length'] / 1609, 1)} mi"         feat["properties"]["tooltip"] = tooltip     return G, nodes, edges 

Уникальное отличие пешеходных маршрутов — в том, что, помимо широты и долготы, при их построении нужно учитывать высоту над уровнем моря. Как и любая хорошая туристическая карта, наша программа должна освоить все топографические особенности местности.

Мой пользовательский API определяет высоты по точным радиолокационным данным 2000 года с шаттла Endeavor. При объёме данных около 18 Гб SRTM GL3 покрывает почти 120 км² суши с шагом в 90 м. Библиотека rasterio позволяет извлекать профили рельефа из отрезков пути, но не загружать в память всю карту:

def get_elevation_profile_of_segment(dataset: rasterio.DatasetReader, coords: list[list]):     """     Get the elevation profile (distance vs. altitude) of a path segment from the list of coordinates.     Args:         dataset: The opened rasterio dataset for the SRTM global topography data.          coords: The path coordinates in [[lon1, lat1], [lon2, lat2], ...] format.      Returns: The distance (in miles) and elevation (in feet) vectors.     """     # coordinates are [lon, lat], flip for rasterio     coords = [[c[1], c[0]] for c in coords]     # convert meters to feet and use rasterio.sample.sample_gen to query each point     elev = [e[0] * 3.28084 for e in sample_gen(dataset, coords)]     d = [0.0]     for j in range(len(coords) - 1):         # use haversine distance         d.append(d[j] + haversine((coords[j][1], coords[j][0]), (coords[j + 1][1], coords[j + 1][0]), Unit.MILES))     return d, elev 

Неплохо, да? Кстати, Google берёт 5 баксов за 1000 запросов к своему Elevation API, но NASA и несколько строчек кода позволяют получить столь же качественные данные совершенно бесплатно.

Dash — моя любимая среда создания интерактивных дашбордов на Python. Её интерфейс состоит из настраиваемых элементов, поэтому при анализе данных можно сосредоточиться на функциях, а не на дизайне. Есть десятки отличных руководств по Dash с примерами, поэтому давайте перейдём сразу к интерфейсу и логике приложения:

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

def layout(place=None):     if place:         global G, G2         G, _, edges = load_map(place)         G2 = nx.Graph(G)          polyline = dl.Polyline(id="route-line", positions=[], color="red", weight=10, fillColor="red", fillOpacity=1.0)         patterns = [dict(repeat='100',                          arrowHead=dict(pixelSize=15, polygon=False, pathOptions=dict(color="red", stroke=True)),                          line=dict(pixelSize=10, pathOptions=dict(color='#f00', weight=20)))]         route_decorator = dl.PolylineDecorator(id="route-arrows", positions=[], patterns=patterns)          url = "https://api.mapbox.com/styles/v1/mapbox/outdoors-v9/tiles/{z}/{x}/{y}?access_token=" + open(".mapbox_token").read()         attribution = '&copy; <a href="https://openstreetmap.org/">OpenStreetMap</a> | &copy; <a href="https://mapbox.com/">Mapbox</a> '          res = html.Div([             html.Div([                 dbc.Label("Choose Mode"),                 dbc.RadioItems(                     options=[                         {"label": "Auto", "value": 1},                         {"label": "Custom", "value": 2}                     ],                     value=1,                     id="mode-select"                 )             ]),             dbc.Label("Distance Target (mi)", html_for="distance-target"),             dbc.Input(id="distance-target", value=10, required=True, type="number", min=0),             dbc.Label("Elevation Minimum (ft)", html_for="elevation-min"),             dbc.Input(id="elevation-min", value=0, required=True, type="number", min=0),             dl.Map([dl.TileLayer(url=url, attribution=attribution),                     polyline,                     dl.GeoJSON(data=edges, zoomToBounds=True, id="trail", hoverStyle=arrow_function(dict(weight=5, color='#666', dashArray=''))),                     # polyline,                     route_decorator                     ],                    id="map",                    style={'width': '100%',                           'height': '40vh',                           'margin': "auto",                           "display": "block",                           "padding-top": "10px"}),             dbc.Button(id="reset-btn", children="Reset", className="btn btn-danger"),             dbc.Button(id="download-btn", children="Export Route"),             dcc.Download(id="download-route"),             # dbc.Button(id="undo-btn", children="Undo"),             html.Br(),             html.P(id="total-distance"),             html.P(id="elevation-gain"),             dcc.Graph(id="profile"),             dcc.Store(id="route-path", data=[]),             dcc.Store(id="map-data", data={"edges": edges})         ],         style={"padding-top": "10px"})         return res     else:         return html.Div("issue") 

После выбора региона сеть троп OSM визуализируется с помощью библиотеки dash-leaflet. При генерации маршрута на двухмерной карте в реальном времени строится и соответствующий профиль высот:

Просмотр трёх направлений петли вокруг горы Олд Рэг в национальном парке Шенандоа (от автора)

Моё приложение позволяет планировать походы в двух режимах. В режиме по умолчанию вы просто нажимаете на пересечения троп на карте и видите, как ваш маршрут «оживает» сегмент за сегментом. В это время программа соединяет точки с помощью алгоритма кратчайшего пути. Вот анимация создания произвольного маршрута по национальному парку Кайахога:

Поэтапное планирование маршрута по национальному парку Кайахога (от автора)

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

• Прирост высоты (постоянно возрастающий).
• Расстояние.
• Движение по как можно более уникальному маршруту (возвращения к пройденным участкам пути на протяжении всего похода сводятся к минимуму).

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

Здесь — это путь, — проходимое расстояние, — целевое расстояние, — дробная доля уникальных участков (если идти по одному пути туда и обратно, то эта доля равна 0,5), — прирост высоты на протяжении всего пути, а — целевой прирост высоты.

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

def find_best_loop(G: nx.Graph, root, target_dist, tol=1.0, min_elev: Optional[Union[int, float]] = None,                    dataset: Optional[rasterio.DatasetReader] = None):     if min_elev and not dataset:         raise ValueError("If asking for elevation data, you must include a rasterio dataset")     error = 1e8     best_path = []     for n in G.nodes():         if nx.has_path(G, root, n):             shortest_path = nx.shortest_path(G, root, n)             paths = nx.all_simple_paths(G, root, n, cutoff=10)             for path in paths:                 if path == shortest_path:                     continue                  path_loop = path + nx.shortest_path(G, n, root, weight="length")[1:]                  path_diversity = len(set(path_loop)) / len(path_loop)                  dist = path_length(G, path_loop)                 e_new = np.abs(dist - target_dist) / path_diversity                 if e_new < error:                     error = e_new                     best_path = path_loop                 if error < tol:                     if min_elev:                         coords_list = path_to_coords_list(G, path_loop)                         _, elev = get_elevation_profile_of_segment(dataset, coords_list)                         elev_gain = elevation_gain(elev)                         if elev_gain < min_elev:                             continue                     return best_path     return best_path 

Клик по карте вблизи пересечения троп определяет корневой узел функции find_best_loop — конечную точку маршрута. Пользователи могут генерировать пешеходные маршруты, которые начинаются и заканчиваются в заданном месте. Это может быть автостоянка или перевалочный пункт. Ниже вы увидите алгоритм, который генерирует суровый марш-бросок с восхождением на Олд Рэг в национальном парке Шенандоа. Здесь выбранные ограничения будут очень кстати.

Используем алгоритм оптимизации и пользовательские ограничения и генерируем идеальный маршрут (от автора)

Ещё одна интересная штука по поводу прироста высоты… Как ни странно, строгого математического определения у него нет. Представляете, что будет, если включать сюда каждый участок с положительным наклоном — даже если он почти плоский? Тогда любая ухабистая дорога на пути существенно повысит это значение. Чтобы сгладить кривые высот, введём вертикальный порог высоты. Это пользовательский гиперпараметр. Он описывается в этой статье.

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

Пользователи также могут экспортировать свои маршруты в GPX (ещё один формат геоданных) и заливать в AllTrails, Gaia, Google Maps и другие приложения с отслеживанием через GPS.

Экспорт маршрута в файл GPX и его загрузка в AllTrails (от автора)

Теперь вы можете выполнять навигацию по вашему личному маршруту с телефона, как по базовому маршруту AllTrails!

Проследил свой поход с помощью мобильного приложения AllTrails (от автора)

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

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

Прокачаться или стать востребованным профессионалом в IT с самого начала помогут наши курсы:

• Профессия Fullstack-разработчик на Python (15 месяцев)
• Профессия Data Scientist (24 месяца)