Sieci przestrzenne

Author

Anna Dmowska, Jakub Nowosad 

library(sf)
library(osmdata)
library(sfnetworks)

1 Podstawowe pojęcia

  • Sieć (Network) jest to zbiór wzajemnie powiązanych obiektów:

    • linii będących krawędziami/segmentami (Edges)
    • punktów będących węzłami (Nodes)

  • Graf

    • wykorzystywany do przedstawiania i badania relacji między obiektami
    • wierzchołki grafu mogą reprezentować jakieś obiekty, natomiast krawędzie mogą wówczas obrazować relacje między takimi obiektami

  • Sieć geoprzestrzenna (ang. geospatial network) to grafy osadzone w przestrzeni geograficznej.

    • wierzchołki i krawędzie w grafie można przedstawić jako obiekty geograficzne
    • wierzchołki reprezentowane są jako punkty o określonych współrzędnych w przestrzeni geograficznej
    • krawędzie mogą być opisane przez indeksy węzłów na ich końcach, np. połączenie lotnicze między dwoma lotniskami (spatially implicit edges)
    • krawędzie mogą być reprezentowane przez obiekt liniowy z przypisaną geometrią, np. sieć dróg (spatially explicit edges)

2 Sieci przestrzenne w R

Dane wykorzystywane w różnego rodzaju analizach powinne być przechowywane w formie uporządkowanej (ang.tidy), tj. tworzyć tabelę, w której:

  • każda zmienna jest przechowywana w osobnej kolumnie,
  • każda obserwacja jest przechowywana w osobnym wierszu,
  • każda komórka zawiera wyłącznie jedną wartość.

Dane sieciowe składające się z węzłów (ang. nodes) połączonych segmentami (ang. edges) nie tworzą uporządkowanego zbioru danych, który łatwo można zapisać w tabeli składającej się z kolumn/wierszy. Można je jednak rozdzielić na dwa obiekty: obiekt przechowujący informacje o węzłach sieci (nodes) oraz obiekt przechowujący informacje o segmentach/krawędziach (edges).

Pakiet sfnetworks konwertuje dane sieciowe do takiej postaci (dwa obiekty przestrzenne, z których jeden opisuje węzły, a drugi krawędzie), a pakiet tidygraph dostarcza narzędzi do przetwarzania tych wirtualnych ramek danych wykorzystując funkcje pakietu dplyr oraz dostarcza narzędzi do obliczania wielu miar grafów.

2.1 Pakiety

  • sfnetworks to pakiet R wykorzystywany do analizy sieci geoprzestrzennych. Pakiet ten łączy funkcjonalności pakietu tidygraph oraz sf. Podstawą strukturą danych tego pakietu jest połączenie dwóch obiektów przestrzennych, z których jeden opisuje węzły (reprezentowane jako punkty), a drugi opisuje krawędzie (reprezentowane jako linie łączące punkty). Tak przechowywane dane stanowią dane wejściowe funkcji dostarczanych zarówno przez pakiet tidygraph jak i sf. Pakiet sfnetworks dostarcza także funkcji dedykowanych do pracy z sieciami geoprzestrzennymi (obliczanie najkrótszej ścieżki, czyszczenie i modyfikacja topologii sieci).

  • tidygraph zapewnia uporządkowany interfejs do przetwarzania danych sieciowych reprezentowanych jako grafy; pozwala m.in na obliczenie wielu miar grafów.

3 Pobieranie danych Open Street Map (OSM)

Jako przykład sieci geoprzestrzennej wykorzystamy dane dotyczące sieci ulic w centrum Poznania. Dane te zostaną pobrane z projektu Open Street Map (OSM), a następnie obiekt wektorowy reprezentujący ulice zostanie przekonwertowany na obiekt reprezentujący sieć geoprzestrzenną używając pakietu sfnetworks.

Dane z projektu OSM mogą być pobrane bezpośrednio do R używając pakietu osmdata. Pakiet osmdata dostarcza narzędzi do pobrania danych dla zdefiniowanego przez użytkownika obszaru. Pobranie danych wymaga:

  • uzyskania informacji o zasięgu obszaru.

Zasięg obszaru może być zdefiniowany za pomocą nazwy (np. Poznan) lub współrzędnych określających zasięg danych (bounding box). Współrzędne muszą być w układzie WGS84. W przypadku, gdy definiujemy zasięg używając nazwy obszaru, musimy wykorzystać funkcję getbb() do uzyskania informacji o zakresie współrzędnych.

getbb("poznan poland")
       min      max
x 16.73159 17.07171
y 52.29192 52.50933
getbb("poznan grunwald poland")
       min      max
x 16.90607 16.90617
y 52.40189 52.40199

Funkcja getbb() pozwala także na zdefiniowanie obszaru w postaci nieregularnego poligonu, a nie prostokątnego zasięgu warstwy

b = getbb ("poznan", format_out = "polygon")
head (b[[1]])
         [,1]     [,2]
[1,] 16.73159 52.46375
[2,] 16.73162 52.46365
[3,] 16.73170 52.46346
[4,] 16.73181 52.46317
[5,] 16.73187 52.46303
[6,] 16.73193 52.46287

Poniższe zapytanie pobierze dane mające atrybut amenity.

q <- opq ("poznan poland") %>%
    add_osm_features(key = "amenity")

Poniższe zapytanie pobierze dane mające atrybut amenity zdefiniowany jako restaurant (w wyniku otrzymamy obiekty typu restauracja).

q <- opq ("poznan poland") %>%
    add_osm_features(key = "amenity", value = "restaurant")

Atrybuty możemy także zdefiniować używając argumentu feautures, w którym definiujemy klucz(key) = wartość(value)

q <- opq ("poznan poland") %>%
    add_osm_features(features = list (
        "amenity" = "university",
        "amenity" = "school"))
  • określenie wynikowego formatu danych używając funkcji osmdata_sf() (np. obiekt klasy sf)
q <- opq ("poznan poland") %>%
    add_osm_features(features = list (
        "amenity" = "university",
        "amenity" = "school")) %>% 
    osmdata_sf()

4 Tworzenie sieci na podstawie OSM

4.1 Zasięg (bounding box)

W poniższym przykładzie zostaną pobrane dane dla fragmentu Poznania wyznaczonego w promieniu 500 m wokół Ratusza. Ratusz ma współrzędne 16.9338 E, 52.4082 N.

Aby wyznaczyć strefę 500 m wokół obiektu punktowego (Ratusza) należy:

  • stworzyć obiekt punktowy klasy sf
p = st_sfc(st_point(st_point(c(16.9338, 52.4082))), crs = 4326)
p
Geometry set for 1 feature 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 16.9338 ymin: 52.4082 xmax: 16.9338 ymax: 52.4082
Geodetic CRS:  WGS 84
POINT (16.9338 52.4082)
  • wyznaczyć strefę bufforową wokół tego punktu
bfr = st_buffer(p, dist = 500)
bfr
Geometry set for 1 feature 
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 16.92632 ymin: 52.40363 xmax: 16.9413 ymax: 52.41275
Geodetic CRS:  WGS 84
POLYGON ((16.92661 52.40926, 16.92656 52.40918,...
  • wyznaczyć zasięg przestrzenny (bounding box) dla strefy bufforowej Zasięg warstwy możemy odczytać używając funkcji st_bbox() z pakietu sf. W wyniku otrzymamy obiekt zawierający 4 współrzędne: xmin, ymin, xmax oraz ymax.
bbox = st_bbox(bfr)
bbox
    xmin     ymin     xmax     ymax 
16.92632 52.40363 16.94130 52.41275

Obiekt ten następnie należy przekształcić na macierz, w której pierwszy wiersz zawiera współrzędne xmin oraz xmax, a drugi wiersz współrzędne ymin oraz ymax. Tak zdefiniowany obiekt możemy użyć w pakiecie osmdata do zdefiniowania zakresu obszaru, dla którego mają zostać pobrane dane.

#convert to matrix
dim(bbox) <- c(2,2)
bbox
         [,1]     [,2]
[1,] 16.92632 16.94130
[2,] 52.40363 52.41275

4.2 Budowa zapytania

Do konstrukcji zapytania mającego na celu pobranie danych z OSM wykorzystuje się funkcje opq(), add_osm_features() oraz osmdata_sf(). W funkcji opq() definiujemy zasięg przestrzenny pobieranych danych, w funkcji add_osm_features() definiujemy atrybuty (features), które mają zostać pobrane, a funkcja osmdata_sf() definiuje wynikowy format danych (dane klasy sf).

W poniższym przykładzie zostaną pobrane dane dla centrum Poznania (zdefiniowanego przez współrzędne w obiekcie bbox) dotyczące sieci dróg dla ruchu pieszego. Pobierzemy drogi (klucz/key highway), typu (wartość/value footway, steps, living street). Dodatkowo foot = yes wskazuje na drogi współdzielone, gdzie część drogi jest przeznaczona dla pieszych.

#Budowa zapytania
library(osmdata)
osm_foot <- opq(bbox) %>% 
  add_osm_features(features = list("highway"="footway",
                                "highway"="steps",
                                "foot"="yes",
                                "highway"="living_street")) %>% 
  osmdata_sf()

4.3 Selekcja obiektów tworzących sieć ulic

Obiekt osm_foot to lista składająca się z kilku typów obiektów: obiektów punktowych, liniowych, poligonowych.

names(osm_foot)
[1] "bbox"              "overpass_call"     "meta"             
[4] "osm_points"        "osm_lines"         "osm_polygons"     
[7] "osm_multilines"    "osm_multipolygons"

Dla celów analizy wybierzemy obiekty liniowe (osm_lines), które następnie możemy zwizualizować.

osm_foot_lines <- osm_foot$osm_lines
osm_foot_lines %>%
  st_geometry() %>% 
  plot()

W następnym kroku zwizualizujemy obiekty zapisane jako poligon (osm_polygons), aby upewnić się, czy w analizowanym obszarze niektóre drogi nie zostały zapisane jako zamknięte linie (zwracane jako poligon).

osm_foot$osm_polygons %>% 
  st_geometry() %>% 
  plot()

Aby uwzględnić powyższe drogi (zapisane jako poligon) w budowanej przez nas sieci ulic musimy najpierw przekonwertować je na obiekty liniowe, a potem dołączyć do wcześniej stworzonego obiektu liniowego (osm_foot_lines).

Funkcja st_cast() z pakietu sf pozwala na konwersję między różnymi typami danych. W przykładzie wykorzystamy tą funkcję do przekonwertowania obiektu poligonowego osm_foot$osm_polygons na typ liniowy (LINESTRING).

Funkcja bind_rows() z pakietu dplyr pozwoli na połączenie dwóch obiektów osm_foot_lines oraz poly_to_lines w jeden.

library(dplyr)
# cast polygons to lines
poly_to_lines <- st_cast(osm_foot$osm_polygons, "LINESTRING")
Warning in st_cast.sf(osm_foot$osm_polygons, "LINESTRING"): repeating
attributes for all sub-geometries for which they may not be constant
# bind all lines together
osm_foot_lines <- bind_rows(osm_foot_lines, poly_to_lines)

Po połączeniu danych otrzymamy sieć ulic, która możemy ponownie zwizualizować.

osm_foot_lines %>% 
  st_geometry() %>% 
  plot()

Do dalszych analiz tabela atrybutów zostanie ograniczona do 4 kolumn: osm_id, name, highway, foot oraz przekonwertowany do PUWG1992 (EPSG: 2180).

osm_foot_lines <- osm_foot_lines %>% 
  select(osm_id, name, highway, foot) %>% 
  st_transform(2180)
write_sf(osm_foot_lines, "data/out_osm_lines.gpkg")

4.4 Wykorzystanie sfnetworks do utworzenia sieci przestrzennej w R

W następnym kroku obiekt liniowy klasy sf osm_foot_lines zostanie przekonwertowany do strkury danych reprezentującej sieć przestrzenną w R. W tym celu wykorzystamy pakiet sfnetworks. Funkcja as_sfnetwork() przekonwertuje obiekt liniowy klasy sf na obiekt klasy sfnetworks reprezentujący sieć przestrzenną. Funkcja ta wymaga podania dwóch argumentów: obiektu wektorowego (osm_foot_lines) oraz zdefiniowania argumentu directed jako FALSE (zostanie utworzony sieć nieskierowana - tj. taka w której krawędzie/segmenty nie mają zdefiniowanego kierunku ruchu).

library(sfnetworks)
foot_net <- as_sfnetwork(osm_foot_lines, directed = FALSE)

Po wyświetleniu obiektu foot_net widzimy, że reprezentuje on teraz sieć składająca się z węzłów (nodes) oraz łączących je krawędzi/segmentów (edges).

plot(foot_net)

Obiekt foot_net jest obiektem klasy sfnetworks. Wpisując nazwę obiektu (foot_net) możemy uzyskać dodatkowe informacje na temat obiektu:

  • liczbę wierzchołków oraz krawędzi w sieci

  • układ współrzędnych

  • typ grafu

    • jest to graf nieskierowany składający się z krawędzi zdefiniowanych przestrznnie (spatially explicit edges) - tj. linii mających przypisaną geometrię.

Obiekt foot_net składa się z dwóch tabel (obiekt typu tibble).

  • Pierwsza tabela przechowuje informacje o węzłach sieci (nodes)

    • kolumna geometria przechowująca współrzędne wierzchołków,

  • Druga tabela przechowuje informacje o krawędziach (edges):

    • atrybut from oraz to to indeks początkowego oraz końcowego węzła, którego współrzędne są zapisane w pierwszej tabeli
    • pozostałe atrybuty zostały pobrane z OSM
foot_net
# A sfnetwork with 1541 nodes and 1130 edges
#
# CRS:  EPSG:2180 
#
# An undirected multigraph with 657 components with spatially explicit edges
#
# A tibble: 1,541 × 1
             geometry
          <POINT [m]>
1 (359112.5 506465.9)
2 (359138.6 506532.7)
3 (359193.2 506680.9)
4   (359274.8 506675)
5 (359023.2 506483.5)
6 (359048.9 506476.2)
# ℹ 1,535 more rows
#
# A tibble: 1,130 × 7
   from    to osm_id   name              highway foot                   geometry
  <int> <int> <chr>    <chr>             <chr>   <chr>          <LINESTRING [m]>
1     1     2 8452538  Aleje Karola Mar… living… <NA>  (359112.5 506465.9, 3591…
2     3     4 8452541  Ignacego Paderew… living… yes   (359193.2 506680.9, 3592…
3     5     6 20910211 Święty Marcin     living… <NA>  (359023.2 506483.5, 3590…
# ℹ 1,127 more rows

W R możemy wskazać, na której tabeli chcemy wykonywać analizy. Służy do tego funkcja activate(), która tymczasowo “wybiera” jedną z tabel, tak aby następnie można było na niej wykonać analizy. W poniższym przykładzie, aktywowaliśmy krawędzie (activate(“edges”)), a następnie używając funkcji mutate() z pakietu dplyr obliczymy długość każdej krawędzi.

foot_net <- foot_net %>% 
  activate("edges") %>% 
  mutate(weight = edge_length())

4.5 Wizualizacja elementów sieci

Aby zwizualizować elementy sieci (węzły lub segmenty) używając pakietu tmap, musimy je przekształcić z powrotem do klasy sf. W tym celu wykorzystujemy funkcję st_as_sf() z pakietu sf, a jako argumenty podajemy obiekt klasy sfnetworks zawierający sieć (w przykładzie foot_net), oraz argument “nodes” jeśli chcemy utworzyć obiekt zawierający same wierzchołki lub “edges” jeśli chcemy utworzyć obiekt zawierający same krawędzie.

library(tmap)
tm_shape(st_as_sf(foot_net, "nodes")) +
  tm_dots(size = 0.3, fill = "blue")

Poniższa mapa przedstawia segmenty/krawędzie, którym przypisano kolor na podstawie atrybutu weight - długości krawędzi.

library(tmap)
tm_shape(st_as_sf(foot_net, "edges")) +
  tm_lines(col = "weight", palette = "YlGn", lwd = 2)

Przestawiając sposób wyświetlania na tmap_mode("view") możemy wyświelić mapę interaktywną. Mapa interaktywna może być przydatna do dokładniejszego przyjrzenia się sieci. tmap_mode("plot") pozwoli na przełączenie sposobu wyświetlania na mapę statyczną.

library(tmap)
tmap_mode("view") # set to interactive mode
tm_shape(st_as_sf(foot_net, "edges")) +
  tm_lines(col = "highway", palette = "Set1", colorNA = "darkgrey", lwd = 2) +
tm_shape(st_as_sf(foot_net, "nodes")) +
  tm_dots(size = 0.3)

4.6 Czyszczenie sieci

Pakiet sfnetworks dostarcza wielu funkcji, które pozwalają na “wyczyszczenie” topologii sieci. Metody te zostały omówione w dokumentacji pakietu https://luukvdmeer.github.io/sfnetworks/articles/sfn02_preprocess_clean.html#network-cleaning-functions

4.6.1 Podział krawędzi

library(tidygraph)
foot_subdivision <- convert(foot_net, to_spatial_subdivision)
Warning: to_spatial_subdivision assumes attributes are constant over geometries

W wyniku tej operacji może zwiększyć się liczba krawędzi.

4.6.2 Usuwanie pseudo-wierzchołków

foot_smooth <- convert(foot_subdivision, to_spatial_smooth)

W efekcie powstanie sieć, która będzie zawierała mniej węzłów.