Questo script scarica, seleziona e plotta i confini delle province italiane estraendo dati da OpenStreetMap

province_italiane.R

library(plyr)
library(osmdata)
library(sf)
library(dplyr)

# Per prima cosa scarico separatamente i confini delle regioni (admin_level = 4), un pezzo alla volta
# Questo mi consente di avere una bounding box più precisa per scaricare successivamente le province
# Se il server da un timeout, bisogna definire più di 4 queries

# Nord Italia (circa)
q1 <- opq(bbox =  c(7.5, 43.7, 11, 46)) %>% 
  add_osm_feature(key = 'admin_level', value = '4') %>% osmdata_sf() %>% 
  unique_osmdata()

# Friuli Venezia Giulia
q2 <- opq(bbox =  c(12, 45, 13, 46.5)) %>% 
  add_osm_feature(key = 'admin_level', value = '4') %>% osmdata_sf() %>% 
  unique_osmdata()

# Centro Italia (circa)
q3 <- opq(bbox =  c(8.2, 40.8, 15.3, 43.2)) %>% 
  add_osm_feature(key = 'admin_level', value = '4') %>% osmdata_sf() %>% 
  unique_osmdata()

# Sud Italia (circa)
q4 <- opq(bbox = c(11.7, 36.6, 18.8, 40.9)) %>% 
  add_osm_feature(key = 'admin_level', value = '4') %>% osmdata_sf() %>% 
  unique_osmdata()

# Metto assieme i risultati delle query
q <- c(q1, q2, q3, q4)
# Cancello i pezzi singoli per liberare un po' di memoria
rm(q1,q2,q3,q4)

# Qui inizia l'ispezione delle query. Ci sono regioni straniere che non c'entrano e vanno tolte.
# Teniamo solo quelle la cui stringa wikipedia contiene 'it:'
regioni_italiane <- q$osm_multipolygons$wikipedia
regioni_italiane <- regioni_italiane[grepl('it:', regioni_italiane)]
regioni <- q$osm_multipolygons %>% filter(wikipedia %in% regioni_italiane)

# Trasformo la tabella in una lista, per poter usare lapply
regioni.list <- split(regioni, seq(nrow(regioni)))

# Poiché la Valle d'Aosta non ha province, va trattata separatamente
# La estraiamo al livello di regione
valle_daosta <- regioni %>% filter(wikipedia == "it:Valle d'Aosta")

# Calcolo la bounding box per ciascuna regione
bbox.list <- lapply(regioni.list, st_bbox)

# Uso bbox.list per fare una lista di query da OSM, questa volta per le province (admin_level = 6)
province <- lapply(bbox.list, function(x) {
  opq(bbox = x) %>% 
    add_osm_feature(key = 'admin_level', value = '6') %>% osmdata_sf() %>% 
    unique_osmdata()
})

# Unisco le query e seleziono i multipolygon
province_unite <- do.call(c, province)
province_mappa <- province_unite$osm_multipolygons

# Ultima selezione, ci sono province straniere anche qua
finale <- province_mappa %>% filter(grepl('it:', wikipedia)) %>%
  filter(!grepl('Distretto di', wikipedia))

# Facciamo un ultimo rbind per reintrodurre la Valle d'Aosta
finale <- rbind.fill(finale, valle_daosta) %>% st_sf()

# Plot finale (statico)
plot(finale['osm_id'])