Zastanawialiście się skąd wziąć listę wszystkich miejsc danego typu – na przykład przedszkola w Warszawie – i narysować je na mapie? Jeśli tak, to z tego postu dowiesz się jak to zrobić. Wpis techniczny, żadna tam analiza. Ale są mapki!
Macie świadomość, że istnieje coś takiego jak Wikipedia, prawda? Zbiór haseł uzupełniany przez ochotników. Istnieje coś podobnego dla map i nazywa się Open Street Map. Ludzie z całego świata dodają informacje o różnych obiektach na mapach – budują drogi, stawiają domki i tak dalej. Dane te można pokazać jako mapy, ale też przeszukiwać. W R mamy pakiet osmdata (przygotowany przez zacną inicjatywę rOpenSci, a jest opakowaniem dla Overpass API), który umożliwi nam skorzystanie z tych informacji. W dalszej części do pokazania (na interaktywnych, przesuwanych i skalowanych) map użyjemy pakietu leaflet (opierający się na bibliotece JavaScript Leaflet, który nota bene również korzysta z Open Street Map.
Ciekawostka na marginesie: ostatnio Google zmienił drastycznie ceny używania Google Maps (ciekawy tekst na ten temat znajdziecie tutaj), zatem korzystanie z Open Street Map może być wybawieniem.
Zaczniemy od wczytania odpowiednich bibliotek:
|
1 2 3 4 5 |
library(tidyverse) # zawsze! głównie dla dplyr i ggplot2 library(sf) # do rysowania danych geograficznych library(osmdata) # clue tego wpisu; używam wersji dev z Githuba, instalacja przez devtools::install_github('osmdatar/osmdata') library(ggmap) # dla statycznych map library(leaflet) # dla map interaktywnych |
W pierwszym kroku zdefiniujmy sobie obszar, z jakiego pobierzemy dane. A że wakacje, to niech to będą Bieszczady:
|
1 |
q0 <- getbb("Wetlina, Poland") %>% opq(timeout = 60) |
Kolejny krok to zdefiniowanie danych, jakie chcemy pobrać.
|
1 |
q1 <- add_osm_feature(q0, key = 'name', key_exact = TRUE) # wszystkie obiekty |
O parametrach funkcji add_osm_feature() będzie za moment, na razie się tym nie przejmujmy. Zobaczmy po prostu co możemy z Open Street Map wyciągnąć, w jakich typach i układach danych.
Na razie zdefiniowaliśmy zapytanie do API, czas pobrać dane. Pobieramy dane w obiektach sf (można użyć również osmdata_sp() jeśli wolicie format sp, ale sf jest wygodniejszy – po szczegóły odsyłam do dokumentacji pakietów sp i sf).
|
1 |
data_sf <- osmdata_sf(q1) |
Dane pobrane (trochę to trwa – w sumie dostajemy jakieś 5GB!). można pobrać też te same dane do lokalnego pliku XML, a później z niego korzystać – odpowiednio modyfikując kod:
|
1 2 |
osmdata_xml(q1, "plik.xml") # pobranie danych do pliku XML data_sf <- osmdata_sf(q1, doc = "plik.xml") # pobranie danych z XMLa |
Co ciekawe – plik na dysku zajmuje mniej niż dane w pamięci – w tym przypadku jest to około 17 MB. Warto ten plik obejrzeć.
Zobaczmy co mamy w pobranych danych:
|
1 |
data_sf |
|
1 2 3 4 5 6 7 8 9 |
Object of class 'osmdata' with: $bbox : 49.0857944,22.4401929,49.1852625,22.5530623 $overpass_call : The call submitted to the overpass API $timestamp : [ Fri 6 Aug 2018 12:36:10 ] $osm_points : 'sf' Simple Features Collection with 188813 points $osm_lines : 'sf' Simple Features Collection with 4316 linestrings $osm_polygons : 'sf' Simple Features Collection with 148 polygons $osm_multilines : 'sf' Simple Features Collection with 18 multilinestrings $osm_multipolygons : 'sf' Simple Features Collection with 30 multipolygons |
Widzimy, że jest to w dużym uproszczeniu lista z kilkoma informacjami, w tym z data.frame osm_points, osm_lines, osm_polygons, osm_multilines czy osm_multipolygons – to są właśnie nasze dane geograficzne.
Na początek narysujmy punkty, ale nie wszystkie 188813 tylko 1% z nich – da nam to pogląd o tym gdzie są te punkty:
|
1 2 3 4 5 |
data_sf$osm_points %>% select(geometry) %>% sample_frac(0.01) %>% ggplot() + geom_sf() |
Pobranie danych w formacie sf pozwala na szybkie ich narysowanie – geom_sf() i po sprawie. Od razu mamy odpowiednie proporcje mapy, odpowiednio podpisane osie itd. Jeśli poszukacie starszych postów i przeanalizujecie kod zobaczycie ile zabawy jest z obiektami sp.
W każdym razie na mapie wyżej widzimy, że mamy nie tylko punkty z okolic Wetliny (jak zdefiniowaliśmy obszar) ale w sumie z całej Polski i nie tylko. To kwestia uporządkowania danych w ramach Open Street Map – dokładnie w to nie wnikałem szczerze mówiąc. Ale możemy przyciąć te dane do tych, które mieszczą się we wskazanym obszarze:
|
1 |
data_sf <- trim_osmdata(data_sf, getbb("Wetlina, Poland", format_out = "polygon")) |
Uwaga – to przycięcie wyrzuca całe obiekty (na przykład linie), które wychodzą poza obszar, a nie przycina ich do obszaru. No i funkcja ma jeszcze problemy – zgłosiłem opdowienie issue do autorów. Po przycięciu mamy:
|
1 |
data_sf |
|
1 2 3 4 5 6 7 8 9 |
Object of class 'osmdata' with: $bbox : 49.0857944,22.4401929,49.1852625,22.5530623 $overpass_call : The call submitted to the overpass API $timestamp : [ Fri 6 Aug 2018 12:36:10 ] $osm_points : 'sf' Simple Features Collection with 5208 points $osm_lines : 'sf' Simple Features Collection with 85 linestrings $osm_polygons : 'sf' Simple Features Collection with 52 polygons $osm_multilines : 'sf' Simple Features Collection with 2 multilinestrings $osm_multipolygons : 'sf' Simple Features Collection with 0 multipolygons |
Widocznie mniej danych, co możemy sprawdzić ponownie rysując punkty (teraz już wszystkie – zrobiło się ich 5208, a więc mniej – było 188813):
|
1 2 3 4 |
data_sf$osm_points %>% select(geometry) %>% ggplot() + geom_sf() |
No dobrze, ale co to za punkty? Gdzie one są fizycznie? Czy da to się zestawić z mapą? Oczywiście. Weźmy mapkę z Google Maps (statyczny obrazek, tutaj taki zawierający zdjęcie satelitarne):
|
1 |
mapa <- get_map(location = "Wetlina, Poland", zoom = 13, maptype = "hybrid") |
i zamiast samych punktów narysujmy polygony, czyli inny typ danych wciągniętych z OSM:
|
1 2 3 |
ggmap(mapa, darken = 0.5) + geom_sf(data = data_sf$osm_polygons, inherit.aes = FALSE, color = "white", size = 1) |
Tak samo możemy narysować linie:
|
1 2 3 |
ggmap(mapa, darken = 0.5) + geom_sf(data = data_sf$osm_lines, inherit.aes = FALSE, aes(color = highway), size = 2) |
czy też te punkty, które mają nazwy:
|
1 2 3 |
ggmap(mapa, darken = 0.7) + geom_sf(data = data_sf$osm_points %>% filter(!is.na(name)), inherit.aes = FALSE, color = "white", size = 0.1) |
Jak widzicie wszystko sprowadza się do tego samego kodu i w zasadzie wybrania odpowiedniego data.frame (elementy listy osm_point, osm_lines itd) i narysowania go poprzez geom_sf().
Ale co to za punkty z nazwami, jakie są te nazwy? Dlaczego wszystkich punktów jest tak dużo, a tych z nazwami mało? Ano dlatego, że każdy polygon czy linia składa się z punktów. Stąd bardzo dużo punktów tworzących linie (na przykład jakieś drogi), a mało punktów konkretnych (typu szkoła, przystanek autobusowy). Jest to dość upierdliwe, ale można się połapać – szczególnie jak poczyta się o samych mapach (OSM oczywiście).
Sprawdźmy punkty z nazwami, może mapka mogłaby być interaktywna? Skorzystajmy tym razem z leaflet. Na początek znajdźmy środek naszego pobranego obszaru, tak aby odpowiednio wycentrować mapkę. Punkty narożne znajdziemy w data_sf$bbox i wystarczy odpowiednio ten ciąg przeliczyć:
|
1 2 3 |
gps_coords <- data_sf$bbox %>% str_split(",") %>% unlist() %>% as.numeric() lng_coord <- (gps_coords[2] + gps_coords[4])/2 lat_coord <- (gps_coords[1] + gps_coords[3])/2 |
Teraz już przygotowujemy mapę:
|
1 2 3 4 5 6 7 |
leaflet() %>% addTiles() %>% setView(lng = lng_coord, lat = lat_coord, zoom = 12) %>% addMarkers(data = data_sf$osm_points %>% filter(!is.na(name)) %>% select(name, geometry), label = ~as.character(name)) |
która wygląda tak:
Przy najechaniu na pinezkę pokazują nam się nazwy punktów – o to nam chodziło. Zwróćcie uwagę, że leaflet też znakomicie radzi sobie z danymi podanymi jako obiekty typu sf (tymi sp też).
Uzbrojeni w podstawowe umiejętności możemy przygotować zestaw funkcji, dzięki którym dalsza eksploracja zawartości Open Street Map będzie sprowadzała się do pojedynczych linii kodu. Zatem po kolei.
Na początek będziemy potrzebowali funkcji, która pobierze nam dane. O parametrach key i value jeszcze będzie, spokojnie:
|
1 2 3 4 5 6 7 8 9 |
getOSMData <- function(f_localisation, f_key, f_value) { data_sf <- getbb(f_localisation) %>% opq() %>% add_osm_feature(key = f_key, value = f_value) %>% osmdata_sf() return(data_sf) } |
Przyda się funkcja, która zwróci nam środek obszaru z pobranymi danymi:
|
1 2 3 4 5 6 7 8 9 10 11 12 |
getLngLat <- function(f_data_sf) { gps_coords <- f_data_sf$bbox %>% str_split(",") %>% unlist() %>% as.numeric() lng_coord <- (gps_coords[2] + gps_coords[4])/2 lat_coord <- (gps_coords[1] + gps_coords[3])/2 return(list(lng = lng_coord, lat = lat_coord)) } |
Teraz przygotujmy komplet funkcji rysujących na mapach interaktywnych (leaflet) odpowiednio punkty, polygony i linie:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 |
plotLeafletMap_points <- function(f_data_sf) { # poniżej używamy kolumny category na mapie - jeśli nie istnieje to ustawiamy, żeby była pusta if(is.null(f_data_sf$osm_points$category)) f_data_sf$osm_points$category <- NA gps_coords <- getLngLat(f_data_sf) leaflet() %>% addTiles() %>% setView(lng = gps_coords$lng, lat = gps_coords$lat, zoom = 11) %>% addMarkers(data = f_data_sf$osm_points, label = ~paste0(if_else(is.na(name), "", as.character(name) ), if_else(is.na(category), "", paste0(" (", category, ")") ) ) ) } plotLeafletMap_polygons <- function(f_data_sf) { gps_coords <- getLngLat(f_data_sf) # workaround konieczny, żeby polygony zadziałały # via https://github.com/ropensci/osmdata/issues/100 names(f_data_sf$osm_polygons$geometry) <- NULL leaflet() %>% addTiles() %>% setView(lng = gps_coords$lng, lat = gps_coords$lat, zoom = 11) %>% addPolygons(data = f_data_sf$osm_polygons, opacity = 1, color = 'red', weight = 1, fillOpacity = 0.2, fillColor ='yellow', smoothFactor = 0.9, label = ~as.character(name)) } plotLeafletMap_lines <- function(f_data_sf) { names(f_data_sf$osm_lines$geometry) <- NULL gps_coords <- getLngLat(f_data_sf) leaflet() %>% addTiles() %>% setView(lng = gps_coords$lng, lat = gps_coords$lat, zoom = 11) %>% addPolylines(data = f_data_sf$osm_lines, color = 'red') } |
Kiedy już mamy gotowy pakiet narzędzi możemy z niego skorzystać.
Na początek sprawdźmy gdzie w Szczecinie są sklepy spożywcze? Potrzebujemy odpowiednich wartości parametrów key i value dla naszej funkcji getOSMData(). Listę znaleźć można na stosownej stronie OpenStreetMap, albo (ale już bez wyjaśnienia co jest czym) korzystając z funkcji available_features() oraz available_tags() z pakietu osmdata.
Niestety jak to w przypadku danych tworzonych przez społeczność musimy liczyć się z tym, że nie są kompletne, mogą nie być aktualne. A tutaj jeszcze dodatkowo – to samo może być opisane różnymi wartościami key i value. Spróbujmy jednak poszukać shop=supermarket czyli supermarketów w Szczecinie:
|
1 |
getOSMData("Szczecin, Poland", f_key = "shop", f_value = "supermarket") %>% plotLeafletMap_points() |
Obiekty uzyskane przez getOSMData() (a precyzyjniej osmdata_sf(), bo nasza funkcja jest tylko opakowaniem dla niej) można łączyć, o na przykład tak (miejsca do spania w Poznaniu):
|
1 2 3 4 5 |
noclegi <- c(getOSMData("Poznań, Poland", 'tourism', 'hotel'), getOSMData("Poznań, Poland", 'tourism', 'motel'), getOSMData("Poznań, Poland", 'tourism', 'hostel'), getOSMData("Poznań, Poland", 'tourism', 'guest_house'), getOSMData("Poznań, Poland", 'tourism', 'apartment')) |
Poszukaliśmy elementów z kategorii (key) tourism w kilku odmianach (hotel, motel, hotel, gościniec i apartament do wynajęcia). Możemy im przypisać teraz wartości w kolumnie category, żeby ładnie na mapie było widać (w dymkach nad pinezkami) co jest czym:
|
1 |
noclegi$osm_points$category <- noclegi$osm_points$tourism |
i na koniec rysujemy sobie całość:
|
1 |
plotLeafletMap_points(noclegi) |
To były punkty, zobaczmy jak działa funkcja rysująca linie – na przykład szukając w Krakowie ścieżek rowerowych:
|
1 2 3 |
c(getOSMData("Kraków, Poland", "route", "bicycle"), getOSMData("Kraków, Poland", "highway", "cycleway")) %>% plotLeafletMap_lines() |
Problem jaki widać od razu to ograniczenie obszaru. O tym między innymi jest w issue 137 do pakietu. W funkcji getOSMData() moglibyśmy dodać odpowiedni fragment z trim_osmdata(), ale to jeszcze nie działa poprawnie. Stąd mamy ścieżkę z Krakowa do Częstochowy czy nawet Wiednia. Poza tym wybraliśmy tutaj tylko dwa typy obiektów wprost nawiązujące do ścieżki rowerowej: route=bicycle oraz highway=cycleway, a kombinacji jest więcej (zobacz tu) i aby było jak należy powinniśmy użyć wszystkich odpowiednich.
Analogicznie do punktów i linii możemy poszukać polygonów. Uniwersytety (właściewie to uczelnie wyższe) w Gdańsku? Proszę bardzo:
|
1 |
getOSMData("Gdańsk, Poland", "amenity", "university") %>% plotLeafletMap_polygons() |
Dane z Open Street Map są stale uzupełniane – zobaczcie na historię ich rozwoju na superowej stronce 10 years of OpenStreetMap. Warto też pooglądać swoją okolicę, może samemu coś dodać? Zerknij na openstreetmap.org, warto. Ciekawy jest również post How to (quickly) enrich a map with natural and anthropic details który wykorzystuje dane pozyskane z pomocą osmdata.
W końcu mapy to świetna zabawa.