honyaki 6 сен 2022 в 21:24

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

10 мин

5.4K

Блог компании SkillfactoryOpen source*Python*Программирование*OpenStreetMap*

Туториал

Перевод

Автор оригинала: Evan Diewald

Как создать приложение с открытым кодом для планирования пеших походов и выбора оптимального маршрута?

Любой заядлый турист знает: всё, что ждёт его в походе, зависит от подготовки. Брать ли дождевики? Сколько идти от одного источника воды до другого? Где лучше ночевать в этой местности в это время года? Но самый важный вопрос звучит проще простого: «Что я там буду делать?» Отвечаем на эти вопросы к старту нашего курса по 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. При генерации маршрута на двухмерной карте в реальном времени строится и соответствующий профиль высот:

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

Движок планирования маршрута

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

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

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

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

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

Формула оптимального пути

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

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

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

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

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

Примечание переводчика о статье по ссылке

В статье по ссылке нет термина «rolling threshold», о котором говорит автор, но есть термин «trackpoint elevation threshold» в разделе «Reducing vertical noise». Дано описание «instead of looking at horizontal distance, it only counts elevation changes once they've passed a certain vertical threshold». Поэтому термин перевёл как «вертикальный порог», так будет просто и понятно

Экспорт маршрута в инструменты навигации

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

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

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

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

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

Пора сворачиваться…

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

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

Ссылки

Aric A. Hagberg, Daniel A. Schult and Pieter J. Swart, «Exploring network structure, dynamics, and function using NetworkX», in Proceedings of the 7th Python in Science Conference (SciPy2008), Gäel Varoquaux, Travis Vaught, and Jarrod Millman (Eds), (Pasadena, CA USA), pp. 11–15, Aug 2008

Map data copyrighted OpenStreetMap contributors and available from https://www.openstreetmap.org

NASA Shuttle Radar Topography Mission (SRTM)(2013). Shuttle Radar Topography Mission (SRTM) Global. Distributed by OpenTopography. https://doi.org/10.5069/G9445JDF. Accessed: 2022–08–30

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

Теги:

Хабы: