OpenStreetMap (OSM) è un progetto collaborativo finalizzato a creare mappe del mondo a contenuto libero. Il progetto punta ad una raccolta mondiale di dati geografici, con scopo principale la creazione di mappe e cartografie.

La caratteristica fondamentale dei dati geografici presenti in OSM è che vengono distribuiti con una licenza libera, la Open Database License: è cioè possibile utilizzarli liberamente per qualsiasi scopo, anche commerciale, con il solo vincolo di citare la fonte e usare la stessa licenza per eventuali lavori derivati dai dati di OSM. (OpenStreetMap)

Per usare i dati di OpenStreetMap, usiamo le funzioni del pacchetto osmdata:

# installazione
install.packages("osmdata")

I pacchetti che servono nel seguito sono:

library(tidyverse)
library(sf)
library(osmdata)

Come prima cosa, cerchiamo i limiti spaziali dei dati che vogliamo utilizzare, la bounding box (bbox), una area definita da due valori di longitudine e latitudine (https://wiki.openstreetmap.org/wiki/Bounding_Box).

# bounding box
bb <- getbb("Teramo")
bb

         min      max
  x 13.34693 14.14632
  y 42.42112 42.89508

I tipi di dati (features) disponibili su OSM sono consultabili alla pagina https://wiki.openstreetmap.org/wiki/IT:Map_Features, oppure attraverso il comando available_features(), e per ciascuna classe, con available_tags().

Se ad esempio siamo interessati ai dati sui luoghi per il tempo libero (feature = leisure), possiamo vedere cosa comprendono con il comando:

available_tags("leisure")

 [1] "adult_gaming_centre" "amusement_arcade"    "bandstand"          
 [4] "beach_resort"        "bird_hide"           "common"             
 [7] "dance"               "disc_golf_course"    "dog_park"           
[10] "escape_game"         "firepit"             "fishing"            
[13] "fitness_centre"      "fitness_station"     "garden"             
[16] "golfcourse"          "hackerspace"         "horse_riding"       
[19] "ice_rink"            "marina"              "miniature_golf"     
[22] "nature_reserve"      "park"                "picnic_table"       
[25] "pitch"               "playground"          "slipway"            
[28] "sports_centre"       "stadium"             "summer_camp"        
[31] "swimming_area"       "swimming_pool"       "track"              
[34] "water_park"         

componiamo la query con la funzione opq() (overpass query), applicata alla bbox del territorio che ci interessa (bb), e aggiungiamo le feature alla quesry (con add_osm_feature("leisure")):

q <- opq(bb) %>% 
  # feature
  add_osm_feature("leisure")

Importiamo i dati in formato sf (vedi Rappresentare i dati geografici):

sf.data <- osmdata_sf(q)

Il risultato è una lista che contiene i seguenti elementi:

names(sf.data)

[1] "bbox"              "overpass_call"     "meta"             
[4] "osm_points"        "osm_lines"         "osm_polygons"     
[7] "osm_multilines"    "osm_multipolygons"

Utilizziamo in questo caso i dati relativi ai punti (sf.data$osm_points), e costruiamo una mappa (vedi voce):

ggplot() +
  geom_sf(data = sf.data$osm_points,
          col = "darkblue") 

I dati che abbiamo scaricato riguardano l'intera area della della bbox, e non i confini del territorio (provincia di Teramo).

Usiamo i confini amministrativi disponibili sul sito dell’Istat (vedi voce):

library(readr)
prov2020 <- st_read("ProvCM01012020_g")

ggplot() +
  # solo la provincia di Teramo, con filter
  geom_sf(data =filter(prov2020, COD_PROV == 67)) +
  geom_sf(data = sf.data$osm_points,
          col = "darkblue") + 
  # tema vuoto
  theme_void()

Per aggiungere delle etichette alla mappa (Vedi voce), possiamo usare uno dei campi del dataset:

names(sf.data$osm_points) %>% head()

  [1] "osm_id"           "name"             "access"           "addr.city"       
  [5] "addr.housenumber" "addr.postcode"   

In questo momento, gli sfondi liberamente scaricabili sono quelli di Stamen. Possono essere scaricati con il pacchetto ggmap, pensato principalmente per le mappe di Google Maps.

library(ggmap)

Scarichiamo lo sfondo stile "terrain", quello tipico delle mappe fisiche, ma senza i nomi ('terrain-background'), utilizzando i limiti dell'area precedente:

ter_map <- get_stamenmap(bb,                # bbox ottenuto con il pacchetto osmdata
              maptype='terrain-background')

Queste mappe possono essere stampate con la funzione ggmap()

ggmap(ter_map)

Per aggiungere questa mappa a quella precedente, faremo prima di tutto attenzione all'ordine (lo sfondo in basso, poi i confini, poi i punti):

ggmap(ter_map) +
  # confini della provincia
  geom_sf(data =filter(prov2020, COD_PROV == 67),
          # area trasparente, perché c'è lo sfondo sotto
          alpha = 0, 
          # colore e dimensione del bordo
          col = "red", size = 1,
          # ignora aes di ggmap
          inherit.aes = FALSE) +

Poiché i sistemi di rappresentazione delle coordinate sono diversi, dobbiamo specificare l'argomento inherit.aes = FALSE per ogni set di dati, e, alla fine, specificare il sistema di riferimento (in questo caso l'EPSG:4326, che sarebbe lo standard, ma non è quello usato da ggmap).

# punti
geom_sf(data = sf.data$osm_points,
        col = "darkblue", size = 0.8,
        inherit.aes = FALSE) + 
  # coordinate reference systems 
  coord_sf(crs=sf::st_crs(4326)) +
  theme_void()

OpenStreetMap.R

# installazione
install.packages("osmdata")
library(tidyverse)
library(sf)
library(osmdata)
 
# bounding box
bb <- getbb("Teramo")
 
# query
q <- # overpass query
   opq(bb) %>% 
   # feature
   add_osm_feature("leisure")
 
# scaricare i dati
sf.data <- osmdata_sf(q)
 
# mappa
ggplot() +
  geom_sf(data = sf.data$osm_points,
          col = "darkblue")
 
# con i dati Istat
library(readr)
prov2020 <- st_read("ProvCM01012020_g")
 
ggplot() +
  # solo la provincia di Teramo, con filter
  geom_sf(data =filter(prov2020, COD_PROV == 67)) +
  geom_sf(data = sf.data$osm_points,
          col = "darkblue") + 
  # tema vuoto
  theme_void()

osm_ggmap.R

library(tidyverse)
library(sf)
library(osmdata)
library(ggmap)
 
# bounding box
bb <- getbb("Teramo")
 
# sfondo
ter_map <- get_stamenmap(bb,               
              maptype='terrain-background')
# mappa
ggmap(ter_map)
 
# mappa e dati OpenStreetMap
ggmap(ter_map) +
  geom_sf(data = sf.data$osm_points,
          col = "darkblue", size = 0.8,
          inherit.aes = FALSE) + 
    # coordinate reference systems 
    coord_sf(crs=sf::st_crs(4326)) +
    theme_void()