Работа с данными OpenStreetMap

Создание картографической основы важно для многих исследований, использующих географические карты городской инфраструктуры. Одной из самых полных баз данных местоположений является некоммерческий проект OpenStreetMap (сокращенно OSM). Здесь мы покажем, как можно самостоятельно загружать карты, отображать на картах объекты, а также изучим вопрос о степени достижимости возможных объектов пожара пожарно-спасательными подразделениями для оценки рисков. Отметим, что особое внимание при оценке расстояния и времени прибытия подразделениями должно уделяться социально значимым объектам и объектам с массовым пребыванием людей. Аналогом библиотек, использованных здесь, является библиотека OSMnx языка Python, которую разрабатывает Geoff Boeing.

Создание картографической основы

Идея создания карты такого рода рассмотрена в блоге, ее автор – Taras Kaduk.

Загрузим необходимые для работы библиотеки. Для получения данных с OpenStreetMap используется библиотека osmdata, для работы с географическими данными на основе стандарта Simple Features – библиотека sf, библиотека osrm является связующим звеном между R и сервисом OSRM для определения расстояния между объектами, времени движения и кратчайшего пути. Следует отметить, что альтернативный подход для доступа к данным, лежащим в основе OpenStreetMap использует библиотека osmdata, а для массовых извлечений больших наборов данных OSM в сжатом pbf-формате – библиотека osmextract.

# общие библиотеки
library(tidyverse)
library(magrittr)

# работа с географическими данными
library(osmdata)
library(osrm)
library(sf)

# шкала масштаба на карте
library(ggspatial)

Зафиксируем географический объект, в нашем случае – город Красноярск.

# здесь можно указать любой другой город
my_place <- "Krasnojarsk Russia"

Сделаем запрос из OSM для нахождения географических границ города.

# загрузим все границы, связанные с объектом
1all_boundaries <- opq(my_place, timeout = 300) |>
2  add_osm_feature(key = "boundary",
                  value = "administrative") |>
3  osmdata_sf() |>
4  unname_osmdata_sf() %>%
  .$osm_multipolygons

1: делает Overpass-запрос; увеличение параметра timeout позволяет делать большие запросы, поскольку время ожидания сервера может истечь до того, как будут доставлены все данные¹;
2: добавляет функции в запрос Overpass, которые задаются парами key-value; подробное описание значений представлено на wiki OSM, для получения полного списка значений можно использовать команду available_features(), в нашем случае это административные границы: областей, районов, округов и т.д.;
3: возвращает sf-объект;
4: удаляет имена из osmdata геометрических объектов.

Выясним, где могут содержаться потенциально возможные границы города.

all_boundaries |> 
  as_tibble() |> 
  dplyr::select(osm_id, name, geom) |> 
  dplyr::filter(grepl("Красноярск", name))

# A tibble: 2 × 3
  osm_id  name                                                              geom
  <chr>   <chr>                                               <MULTIPOLYGON [°]>
1 190090  Красноярский край          (((111.3179 73.85324, 111.3222 73.85231, 1…
2 1157393 городской округ Красноярск (((93.01595 56.06642, 93.01739 56.06589, 9…

Выделим из полученного множества границ только нужную нам, соответствующую границе городского округа.

boundary_KRSK <- all_boundaries |> 
                   dplyr::filter(osm_id == 1157393)

Другой простой способ получения границы города – воспользоваться библиотекой nominatimlite для (обратного) геокодирования на основе Nominatim API.

boundary_KRSK <- nominatimlite::geo_lite_sf(my_place, points_only = FALSE)

Теперь получим картографические данные, включающие в себя дорожную сеть улиц, железные дороги и водные объекты – реки и озера. Фиксируем значения ключа highway.

Код: типы загружаемых дорог

# типы загружаемых дорог и толщина линий на карте
highway_sizes <- tibble::tribble(
  ~highway, ~highway_group,
  "motorway",      "large",
  "trunk",         "large",
  "motorway_link", "large",
  "primary",       "large",
  "primary_link",  "large",
  "secondary",     "medium",
  "secondary_link","medium",
  "tertiary",      "medium",
  "tertiary_link", "medium",
  "trunk_link",    "medium",
  "residential",   "small",
  "living_street", "small",
  "unclassified",  "small",
  "service",       "small",
  "road",          "small",
)

Загрузим дорожную сеть города.²

# streets (улицы)
streets_osm <- opq(my_place, timeout = 300) |>
  add_osm_feature(key = "highway", 
                  value = highway_sizes$highway) |>
  osmdata_sf() |> 
  unname_osmdata_sf()

Выделим только типы линий osm_lines и оставим только улицы, находящиеся внутри границы boundary_KRSK.

streets <- streets_osm$osm_lines %>%
  mutate(length = as.numeric(st_length(.))) |>
  left_join(highway_sizes, by = "highway") |> 
  st_intersection(boundary_KRSK)

streets |>
  dplyr::select(osm_id, name, highway, 
                maxspeed, oneway, surface) |>
  head(5)

Simple feature collection with 5 features and 6 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: 92.9288 ymin: 56.0738 xmax: 92.94057 ymax: 56.08547
Geodetic CRS:  WGS 84
    osm_id             name highway maxspeed oneway surface
1 17639885   Северное шоссе primary       60    yes asphalt
2 25375023 Енисейский тракт primary       60    yes asphalt
3 25375024      улица 9 Мая primary       60    yes asphalt
4 25375027 Енисейский тракт primary       60    yes asphalt
5 25375030 Енисейский тракт primary       60    yes asphalt
                            geom
1 MULTILINESTRING ((92.93389 ...
2 MULTILINESTRING ((92.93911 ...
3 MULTILINESTRING ((92.93978 ...
4 MULTILINESTRING ((92.93999 ...
5 MULTILINESTRING ((92.93936 ...

Загрузим данные соответствующие железным дорогам.

Код: железные дороги

# railways (железные дороги)
railways_osm <- opq(my_place, timeout = 300) |>
  add_osm_feature(key = "railway", 
                  value = "rail") |>
  osmdata_sf() |> 
  unname_osmdata_sf()

railways <- railways_osm$osm_lines |> 
  dplyr::select()

Получим данные по водным объектам.

Код: водные объекты

# water (водные объекты)
water_osm <- opq(my_place, timeout = 300) |>
  add_osm_feature(key = "natural", 
                  value = "water") |>
  osmdata_sf() |> 
  unname_osmdata_sf()

river_osm <- opq(my_place, timeout = 300) |>
  add_osm_feature(key = "waterway", 
                  value = c("river", 
                            "riverbank")) |>
  osmdata_sf() |> 
  unname_osmdata_sf()

water <- 
  c(water_osm, river_osm) %>%
  .$osm_multipolygons |>
  dplyr::select(osm_id, name)

Создадим базовую карту, нанеся на нее соответствующие объекты.

Код: базовая карта

base_map <- ggplot() +
  # дорожная сеть
  geom_sf(data = streets |>
          dplyr::filter(highway_group == "large"),
          linewidth = 0.4,
          color = "grey30") +
  geom_sf(data = streets |>
          dplyr::filter(highway_group == "medium"),
          linewidth = 0.2,
          color = "grey35") +
  geom_sf(data = streets |>
          dplyr::filter(highway_group == "small"),
          linewidth = 0.1,
          color = "grey40") +
  # железные дороги
  geom_sf(data = railways,
          color = "grey30",
          linewidth = 0.3,
          linetype = "dotdash",
          alpha = 0.6) +
  # водные объекты
  geom_sf(data = water,
          fill = "steelblue",
          lwd = 0,
          alpha = 0.3) +
  hrbrthemes::theme_ipsum() +
  theme(plot.margin = ggplot2::margin(0.01, 0.01, 0.01, 0.01, "cm"))

Карта-подложка готова, теперь можно использовать ее элементы для отображения всего города или его части. Варьируя границы, мы можем получить карту достаточно подробного масштаба.

base_map +
  # границы города
  geom_sf(data = boundary_KRSK, 
          color = "red",
          alpha = 0.1) +
  # ограничения
  coord_sf(xlim = c(92.64, 93.19), 
           ylim = c(55.90, 56.14),
           expand = FALSE)

Дополнительные данные из OpenStreetMap

Напомним, что базовая карта состоит из трех основных составляющих: streets, railways и water. Мы можем использовать информацию из OpenStreetMap для дальнейшего анализа. В частности, переменная speed объекта streets отвечает за максимальную разрешенную скорость для различных транспортных средств.

Код: максимальные разрешенные скорости движения

# максимальные разрешенные скорости движения из OpenStreetMap
streets_speed <- streets |> 
  # рассмотрим только данные без пропусков
  dplyr::filter(!is.na(maxspeed)) |> 
  # заменим коды значениями из OSM
  mutate(maxspeed_new = 
           case_when(
             maxspeed == "RU:living_street" ~ "20",
             maxspeed == "RU:rural" ~ "90",
             maxspeed == "RU:urban" ~ "60",
             .default = maxspeed
             )
           ) 

# преобразуем переменные
streets_speed$maxspeed_new <- 
  as.numeric(streets_speed$maxspeed_new)

# создадим факторную переменную для скорости
streets_speed$max_speed_factor <-
  factor(streets_speed$maxspeed_new,
         levels = c(5, 8, 10, 20, 30, 40, 49, 50, 60, 70, 90),
         labels = c("5 км/ч", "8 км/ч", 
                    "10 км/ч", "20 км/ч", 
                    "30 км/ч", "40 км/ч", 
                    "49 км/ч", "50 км/ч", 
                    "60 км/ч", "70 км/ч", 
                    "90 км/ч"))

streets_speed |> 
  st_drop_geometry() %>% 
  count(., max_speed_factor, sort = TRUE)

   max_speed_factor    n
1           60 км/ч 2556
2           20 км/ч  394
3           40 км/ч  357
4            5 км/ч   83
5           30 км/ч   26
6           90 км/ч   22
7           10 км/ч    6
8           49 км/ч    1
9           50 км/ч    1
10          70 км/ч    1
11           8 км/ч    1

Код: максимальные разрешенные скорости движения (карта)

# карта: максимальная разрешенная скорость движения
base_map +
  # дорожная сеть
  geom_sf(data = streets_speed,
          linewidth = 0.67,
          aes(color = max_speed_factor)) +
  viridis::scale_color_viridis(discrete = TRUE,
                               option = "turbo") +
  # географические границы части города
  coord_sf(xlim = c(92.8002, 93.0407), 
           ylim = c(55.9823, 56.1003),
           expand = FALSE) +
  labs(color = "максимальная \nразрешенная \nскорость") +
  theme_void() +
  theme(text = element_text(size = 14))

Как видно из рисунка, база OpenStreetMap содержит неполную информацию по каждому участку дорог. Аналогично можно рассмотреть и визуализировать другие характеристики: дороги с односторонним движением (переменная oneway), тип дорожного покрытия (переменная surface), количество полос движения (переменная lanes) и т.д.

Нанесение на карту города различных объектов и ограниченных зон

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

# получение данных по районам города из OSM
districts_osm <- opq(my_place) |>
  add_osm_feature(key = "admin_level", value = 9) |>
  osmdata_sf() |>
  unname_osmdata_sf()

Получим данные по районам, также выделив площадь каждого района.³

# границы районов с площадями
districts <- districts_osm %>% 
  .$osm_multipolygons |>
  dplyr::select(osm_id, name) %>% 
  mutate(area = st_area(.),
         region_area = as.numeric(area))

Найдем центроиды районов.

# центроиды районов
regions <- cbind(districts, 
                 st_coordinates(st_centroid(districts))) |>
  rename(c_lon = X, c_lat = Y)

# районы города
regions

Simple feature collection with 7 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 92.62752 ymin: 55.91184 xmax: 93.16844 ymax: 56.13382
Geodetic CRS:  WGS 84
   osm_id                  name     area region_area    c_lon    c_lat
1 5854088 Железнодорожный район 11472481    11472481 92.82686 56.01383
2 5854089       Кировский район 23189104    23189104 92.96615 55.98596
3 5854090       Ленинский район 52892774    52892774 93.02779 56.02276
4 5854091     Октябрьский район 90609016    90609016 92.73413 56.01672
5 5854092    Свердловский район 75824721    75824721 92.84575 55.96372
6 5854093       Советский район 88194468    88194468 92.97927 56.07533
7 5854094     Центральный район 34350926    34350926 92.88060 56.04314
                            geom
1 MULTIPOLYGON (((92.84876 56...
2 MULTIPOLYGON (((92.96183 55...
3 MULTIPOLYGON (((93.04703 56...
4 MULTIPOLYGON (((92.66951 55...
5 MULTIPOLYGON (((92.96183 55...
6 MULTIPOLYGON (((93.09394 56...
7 MULTIPOLYGON (((92.84994 56...

Приведем карту районов города.⁴

Код: карта районов города

# удалим "район" из названия
regions$name <- str_remove(regions$name, " район")

set.seed(2024)

base_map +
  # районы города
  geom_sf(data = regions, 
          aes(fill = region_area / 1000000), 
          color = "black", 
          alpha = 0.5) + 
  # цвет
  viridis::scale_fill_viridis(option = "plasma", 
                              direction = 1) +
  labs(fill = "площадь района (в кв. км):", 
       x = "", y = "") +
  theme(legend.key.size = unit(0.4,"cm"),
        legend.key.width = unit(1.8,"cm"),
        legend.position = "top") +
  # названия районов
  ggrepel::geom_text_repel(data = regions, 
                           aes(c_lon, c_lat,
                               label = name, 
                               color = I("black")), 
                           color = "white",     
                           bg.color = "grey30", 
                           bg.r = 0.15,          
                           size = 5, alpha = 0.9) +
  # ограничения
  coord_sf(xlim = c(92.64, 93.19), 
           ylim = c(55.90, 56.14),
           expand = FALSE)

Аналогично можно загрузить и отобразить на карте информацию по зданиям и сооружениям.

# buildings (здания)
buildings_osm <- opq(my_place, timeout = 300) |>
  add_osm_feature(key = "building") |>
  osmdata_sf()

buildings_polygons <- buildings_osm %>% 
  .$osm_polygons

buildings_multipolygons <- buildings_osm %>% 
  .$osm_multipolygons

buildings <- bind_rows(buildings_polygons, 
                       buildings_multipolygons)

Рассмотрим количество объектов, которые принадлежат улицам города.

count(buildings |> 
        st_drop_geometry() |> 
        as_tibble(), 
      addr.street, sort = TRUE) |>
  na.omit()

# A tibble: 1,016 × 2
   addr.street                                  n
   <chr>                                    <int>
 1 проспект им. газеты Красноярский Рабочий   457
 2 Центральная улица                          401
 3 Советская улица                            337
 4 Лесная улица                               334
 5 улица Калинина                             310
 6 улица 60 лет Октября                       232
 7 улица Карла Маркса                         231
 8 Новая улица                                221
 9 Базайская улица                            220
10 улица Березина                             218
# ℹ 1,006 more rows

Код: улицы города

base_map +
  geom_sf(data = buildings, linewidth = 0.15, 
          aes(fill = `addr:street`, 
              color = `addr:street`)) +
  # географические границы части города
  coord_sf(xlim = c(92.8492, 93.0487), 
           ylim = c(55.9863, 56.0603),
           expand = FALSE) +
  theme_void() +
  theme(legend.position = "none")

Другой пример – построим карту, учитывающую этажность зданий.

Код: этажность зданий

# этажность
base_map +
  geom_sf(data = buildings |>
            dplyr::filter(as.numeric(`building:levels`) >= 0) |>
            mutate(`building:levels` = factor(`building:levels`, 
                                              levels = c(0:26, 29, 30))), 
          linewidth = 0.1,
          aes(fill = `building:levels`,
              color = `building:levels`)) +
  # географические границы города
  coord_sf(xlim = c(92.88, 92.95), 
           ylim = c(56.02, 56.06),
           expand = FALSE) +
  # шкала аннотаций
  annotation_scale(location = "br", 
                   width_hint = 0.5, 
                   style = "ticks") +
  theme(legend.position = "right") +
  viridis::scale_fill_viridis(option = "turbo", 
                              discrete = T) +
  labs(fill = "этажность") +
  guides(color = "none") +
  hrbrthemes::theme_ipsum()

Анализ зон транспортной доступности

Рассмотрим транспортную достижимость объектов по времени движения по транспортной сети. Задачи, связанные с построением зон доступностей, рассмотрены, например, в главе Сетевой анализ курса [Т.Е. Самсонов]. В частности, Федеральный закон от 22.07.2008 № 123-ФЗ (ред. от 25.12.2023) “Технический регламент о требованиях пожарной безопасности” в статье 76 Требования пожарной безопасности по размещению зданий пожарных депо на территориях населенных пунктов говорит о том, что здания пожарных депо на территориях населенных пунктов следует размещать исходя из условия, что время прибытия первого подразделения к месту вызова в городских населенных пунктах не должно превышать 10 минут, в сельских населенных пунктах 20 минут. В связи с этим, актуальным является вопрос достижимости объектов с массовым пребыванием людей (т.е. зданий или сооружений (кроме жилых домов), в котором может одновременно находиться 50 и более человек) из пожарно-спасательных подразделений.

Точки исследования

Помимо карты-подложки нам понадобятся объекты исследования, пусть это будут пожарно-спасательные подразделения, дислоцируемые на территории местного пожарно-спасательного гарнизона (сокращенно ПСЧ). Сформируем точки, которые им соответствуют, это можно сделать случайным образом, например, с помощью команды sf::st_sample(), либо приписать заданным значениям как в нашем случае – центроидам районов города.

# пожарные части
# точки выбраны как центроиды районов
fire_stations <-
regions |>
  mutate(name = c(paste("ПСЧ", 1:7, sep = "-"))) |>
  dplyr::select(name, c_lon, c_lat) |>
  st_drop_geometry() |>
  as_tibble()

fire_stations

# A tibble: 7 × 3
  name  c_lon c_lat
  <chr> <dbl> <dbl>
1 ПСЧ-1  92.8  56.0
2 ПСЧ-2  93.0  56.0
3 ПСЧ-3  93.0  56.0
4 ПСЧ-4  92.7  56.0
5 ПСЧ-5  92.8  56.0
6 ПСЧ-6  93.0  56.1
7 ПСЧ-7  92.9  56.0

Отметим, что в картографии объекты должны рассматриваться в соответствующей географической проекции, которая учитывает искажения Земли. В R библиотека sf позволяет сделать перевод в географическую систему координат CRS (coordinate reference system).

# CRS (coordinate reference system) проекция в EPSG:4326
fire_stations_sf <- fire_stations |>
  st_as_sf(coords = c("c_lon", "c_lat"), 
           crs = 4326)

fire_stations_sf

Simple feature collection with 7 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 92.73413 ymin: 55.96372 xmax: 93.02779 ymax: 56.07533
Geodetic CRS:  WGS 84
# A tibble: 7 × 2
  name             geometry
* <chr>         <POINT [°]>
1 ПСЧ-1 (92.82686 56.01383)
2 ПСЧ-2 (92.96615 55.98596)
3 ПСЧ-3 (93.02779 56.02276)
4 ПСЧ-4 (92.73413 56.01672)
5 ПСЧ-5 (92.84575 55.96372)
6 ПСЧ-6 (92.97927 56.07533)
7 ПСЧ-7  (92.8806 56.04314)

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

Код: диаграммы Вороного

# диаграммы Вороного
base_map +
  geom_sf(data = fire_stations_sf |>
            st_combine() |>
            st_voronoi() |> 
            st_collection_extract() |> 
            st_intersection(boundary_KRSK), 
          color = "darkblue",
          linewidth = 0.7,
          alpha = 0.01) +
  # точки с ПСЧ
  geom_sf(data = fire_stations_sf,
          color = "magenta4", 
          alpha = 1, 
          shape = 8, 
          size = 3, 
          stroke = 1.2) +
  # маркеры
  ggrepel::geom_label_repel(data = fire_stations, 
                            aes(c_lon, c_lat, label = name), 
                            fontface = "bold",
                            size = 4, alpha = 0.9) +
  # географические границы части города
  coord_sf(xlim = c(92.72, 93.04), 
           ylim = c(55.96, 56.09),
           expand = FALSE) +
  # шкала аннотаций
  annotation_scale(location = "tl", 
                   width_hint = 0.5, 
                   style = "ticks") +
  labs(x = "", y = "")

Выберем учебные заведений города Красноярска в качестве примера объектов исследования и проанализируем, все ли из этих зданий находятся в 10-минутной зоне доступности из ПСЧ.

# школы, находящиеся внутри границы города
schools <- buildings |> 
  dplyr::filter(building %in% c("school")) |> 
  st_intersection(boundary_KRSK)

Вместо самих зданий мы выберем их центроиды.

# центроиды объектов
schools_points_KRSK <- st_coordinates(st_centroid(schools)) |>
  rename(sch_lon = X, sch_lat = Y) |>
  dplyr::select(sch_lon, sch_lat) |>
  st_drop_geometry() |>
  st_as_sf(coords = c("sch_lon", "sch_lat"), crs = 4326)

Достижимость объектов на основе матрицы расстояний

Получим матрицу расстояний между ПСЧ и учебными заведениями города с использованием библиотеки osrm, являющейся интерфейсом между R и сервисом маршрутизации OSRM на основе OSRM API.

# матрица расстояний между объектами
distancetable = osrmTable(src = fire_stations_sf, 
                          dst = schools_points_KRSK)

Преобразуем таблицу расстояний.

schools_distances <-
  schools_points_KRSK |>
  mutate(mintime = distancetable$durations |>
           as_tibble() |> 
           summarise(across(where(is.numeric), min)) |>
           t() %>% 
           as.vector(.))

Отметим, что альтернативным подходом для нахождения матрицы, которая отражает расстояния либо время движения между объектами является использование библиотеки r5r – интерфейсу к сервису $R^5$, что расшифровывается как Rapid Realistic Routing on Real-world and Reimagined networks и библиотеки accessibility, которые взяты в основу книги [R.H.M. Pereira, D. Herszenhut].

Нанесем на карту результат, выделив цветом расчетное время прибытия и формой значка те учебные заведения, где время прибытия превысит расчетные 10 минут (треугольник).

Код: достижимость объектов

set.seed(2024)
base_map + 
  # нанесение школ
  geom_sf(data = schools_distances, 
          aes(fill = mintime, 
              # форма значка
              shape = I(ifelse(mintime >= 10, 24, 22))),
          alpha = 0.95, size = 3.7
  ) +
  # цвет
  viridis::scale_fill_viridis(option = "turbo", 
                              breaks = c(1, 5, 10, 15),
                              labels = c("1 мин", "5 мин", 
                                         "10 мин", "15 мин"),
           limits = c(0, max(schools_distances$mintime))) +
  # точки с ПСЧ
  geom_sf(data = fire_stations_sf,
          color = "magenta4", 
          alpha = 1, 
          shape = 8, 
          size = 3, 
          stroke = 1.2)  +
  # маркеры
  ggrepel::geom_label_repel(data = fire_stations, 
                            aes(c_lon, c_lat, label = name), 
                            fontface = "bold",
                            size = 4, alpha = 0.9) +
  # географические границы части города
  coord_sf(xlim = c(92.79, 93.04), 
           ylim = c(55.97, 56.08),
           expand = FALSE) +
  # шкала аннотаций
  annotation_scale(location = "tl", 
                   width_hint = 0.5, 
                   style = "ticks")  +
  # оформление темы
  theme(legend.key.size = unit(0.5,"cm"),
        legend.key.width = unit(2,"cm"),
        legend.position = "top",
        plot.margin = ggplot2::margin(0.01, 0.01, 0.01, 0.01, "cm")) +
  labs(x = "", y = "",
       fill = "время прибытия:")

Очевидно, что не все учебные заведения города при данной условной расстановке ПСЧ являются достижимыми за отведенное время 10 минут.

schools_distances |>
  mutate(
    late_category = case_when(
      mintime >  10 ~ "> 10 мин",
      mintime <= 10 ~ "<= 10 мин"
    )
  ) |>
  janitor::tabyl(late_category) |>
  janitor::adorn_pct_formatting(digits = 0) |>
  purrr::set_names("категория", "количество", "процент")

 категория количество процент
 <= 10 мин        154     81%
  > 10 мин         36     19%

Несмотря на то, что мы решили задачу по анализу и визуализации доступности объектов относительно выбраанных объектов, наше решение обладает рядом существенных недостатков. Во-первых, для большого количества объектов мы не сможем осуществить большой объем запросов (максимум – один запрос в секунду, никакого веб-скреппинга и т.д.). Кроме того, здесь мы не можем самостоятельно настраивать скорости движения автомобиля по различным типам дорог.

Одним из оптимальных решений является упрощение дорожной сети с помощью библиотек sfnetworks, dodgr, либо OSMnx в Python до графа дорожной сети, что делает настройки более гибкими и позволяет использовтаь алгоритмы теории графов. Этот подход будет рассмотрен в другой статье.

Построение оптимальных маршрутов

Покажем, как можно строить оптимальные маршруты с помощью библиотеке osrm. Выберем две произвольных точки, которые будут соответствовать учебному заведению и пожарной части.

set.seed(2024)
fire_stations_sample <- fire_stations_sf |>
  sample_n(1)

schools_sample <- schools_points_KRSK |>
  mutate(id = row_number()) |>
  sample_n(1)

Построим оптимальный маршрут движения между от первой до второй выбранной точки пользуясь командой osrmRoute.

route <- osrmRoute(src = fire_stations_sample, 
                   dst = schools_sample,
                   overview = "full")

route

Simple feature collection with 1 feature and 4 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 92.78692 ymin: 55.97694 xmax: 92.96585 ymax: 56.01863
Geodetic CRS:  WGS 84
  src dst duration distance                           geom
1   1   1    22.28  16.3655 LINESTRING (92.96585 55.985...

Мы видим, что команда osrmRoute() нашла расчетное время движения duration и длинуdistance маршрута. На карте ниже можно увидеть построенный маршрут.

Код: построение маршрута

# построение маршрута
set.seed(2024)
base_map + 
  # здания
  geom_sf(data = buildings, linewidth = 0.05) +
  # маршрут
  geom_sf(data = route, 
          color = "red", 
          linewidth = 1) +
  # ПСЧ -- исходная точка
  geom_sf(data = fire_stations_sample,
          color = "magenta4", 
          alpha = 1, 
          shape = 8, 
          size = 3, 
          stroke = 1.3) +
  # школа -- конечная точка
  geom_sf(data = schools_sample,
          color = "darkblue", 
          alpha = 1, 
          shape = 8, 
          size = 3,
          stroke = 1.3) +
  # ярлык
  geom_sf_label(data = fire_stations_sample, 
                aes(label = name),
                nudge_y = -0.003, 
                fontface = "bold",
                size = 4, alpha = 0.9) +
  # географические границы города
  coord_sf(xlim = c(92.77, 92.99), 
           ylim = c(55.97, 56.03),
           expand = FALSE) +
  theme_void() +
  labs(x = "", y = "")

Заметим, что фукция osrmIsochrone() в osrm позволяет находить изохроны (области, ограниченные фиксированным временем движения), тем не менее, она обладает теми же недостатками, о которых говорилось выше; кроме того, область построения изохроны может быть неодносвязной, некорректно вычисляться и т.д., даже несмотря на варьирования параметра res отвечающего за разрешение.

Заключение

В статье были кратко рассмотрены следующие вопросы в применении к географическим данным городской инфраструктуры:

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

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

Существует множество инструментов для решения такого рода задач как на языке R, так и на языке Python, и библиотека osrm является не единственной библиотекой в этом классе.

Session Info

─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.3.2 (2023-10-31)
 os       macOS Monterey 12.1
 system   aarch64, darwin20
 ui       X11
 language (EN)
 collate  ru_RU.UTF-8
 ctype    ru_RU.UTF-8
 tz       Asia/Krasnoyarsk
 date     2024-03-01
 pandoc   2.18 @ /Users/materov/opt/miniconda3/envs/ox/bin/ (via rmarkdown)
 quarto   1.4.550

─ Packages ───────────────────────────────────────────────────────────────────
 package     * version date (UTC) lib source
 dplyr       * 1.1.4   2023-11-17 [1] CRAN (R 4.3.1)
 forcats     * 1.0.0   2023-01-29 [1] CRAN (R 4.3.0)
 ggplot2     * 3.5.0   2024-02-23 [1] CRAN (R 4.3.1)
 ggspatial   * 1.1.9   2023-08-17 [1] CRAN (R 4.3.0)
 lubridate   * 1.9.3   2023-09-27 [1] CRAN (R 4.3.1)
 magrittr    * 2.0.3   2022-03-30 [1] CRAN (R 4.3.0)
 osmdata     * 0.2.5   2023-08-14 [1] CRAN (R 4.3.0)
 osrm        * 4.1.1   2023-03-30 [1] CRAN (R 4.3.0)
 purrr       * 1.0.2   2023-08-10 [1] CRAN (R 4.3.0)
 readr       * 2.1.5   2024-01-10 [1] CRAN (R 4.3.1)
 sessioninfo * 1.2.2   2021-12-06 [1] CRAN (R 4.3.0)
 sf          * 1.0-15  2023-12-18 [1] CRAN (R 4.3.1)
 stringr     * 1.5.1   2023-11-14 [1] CRAN (R 4.3.1)
 tibble      * 3.2.1   2023-03-20 [1] CRAN (R 4.3.0)
 tidyr       * 1.3.1   2024-01-24 [1] CRAN (R 4.3.1)
 tidyverse   * 2.0.0   2023-02-22 [1] CRAN (R 4.3.0)

 [1] /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library

──────────────────────────────────────────────────────────────────────────────

Сноски

Можно также установить глобальную опцию options(timeout = max(300, getOption("timeout"))).↩︎
Данные по дорожной сети из OSM могут содержать некорректные объекты, например, объектам могут соответствовать неверные координаты, что может приводит к смещениям объектов. В этом случае все неточности можно исправить с помощью функций библиотеки sf.↩︎
Отметим, что для некоторых населенных пунктов также следует рассмотреть .$osm_polygons при нахождении границ районов.↩︎
Для некоторых городов (например, Новосибирска) районы могут быть разделены на несколько компонентов, что нужно учитывать при нахождении центроидов и нанесение их на карты.↩︎

Ссылка для цитирования

BibTeX

@online{матеров2021,
  author = {Матеров, Е.Н.},
  title = {Анализ географически распределенных объектов городской
    инфраструктуры},
  date = {2021-11-22},
  url = {https://www.naukaidannye.netlify.app/blog/posts/2021-11-22-spatial},
  langid = {ru}
}

На публикацию можно сослаться как

Матеров Е.Н. Анализ географически распределенных объектов городской инфраструктуры [Electronic resource]. 2021. URL: https://www.naukaidannye.netlify.app/blog/posts/2021-11-22-spatial.