In questo post vedremo come creare una choropleth map, ossia una mappa geografica con le diverse regioni colorate in base ad un certo valore numerico.
Per far questo, ecco ciò di cui abbiamo bisogno:
** libreria
sp** libreria
Hmisc** la mappa bianca: su questo sito è possibile ottenere le mappe in un formato leggibile da R
** i dati da plottare: per questo esempio i dati sono stati presi dal sito Istat
Mettiamoci al lavoro.
Innanzitutto vediamo come caricare in R la mappa bianca, da colorare successivamente. Identifichiamo la nazione di nostro interesse sul sito Gadm, copiamo il link e importiamolo in R con questo codice:
con <- url("http://gadm.org/data/rda/ITA_adm1.RData")
print(load(con))
close(con)
Il primo problema che si pone è quello di capire in che modo verranno colorate le regioni. Con la funzione
spplot della library(sp), sarà sufficiente creare un vettore di colori, ciascuno dei quali verrà assegnato alla singola regione. In che ordine vengono presentate le regioni? Richiamiamolo con il seguente codice:
summary(gadm$NAME_1)
Abruzzo Apulia Basilicata
1 1 1
Calabria Campania Emilia-Romagna
1 1 1
Friuli-Venezia Giulia Lazio Liguria
1 1 1
Lombardia Marche Molise
1 1 1
Piemonte Sardegna Sicily
1 1 1
Toscana Trentino-Alto Adige Umbria
1 1 1
Valle d'Aosta Veneto
1 1
Osservate che le regioni sono state messe in ordine alfabetico, ad eccezione della Puglia che viene chiamata Apulia. Quando assegneremo i colori alle regioni, dobbiamo seguire questo ordine.
A questo punto dal sito Istat cerchiamo i valori da analizzare. In questo esempio, andremo a dividere le regioni italiana in base alla percentuale di donne nubili. Cerchiamo i valori sul sito e importiamoli in R (clickate su view plain e allargata la finestra per visualizzare correttamente questa tabella):
data <- read.table("clipboard", header=T)
data
Regioni Maschi_Celibi Totale_Maschi Femmine_Nubili Totale_Femmine Maschi_Femmine
1 Abruzzo 284385 648680 242027 685995 1334675
2 Puglia 883027 1979254 809566 2100448 4079702
3 Basilicata 130276 289275 112162 301326 590601
4 Calabria 452617 978789 398402 1029920 2008709
5 Campania 1324589 2820078 1207237 2992884 5812962
6 Emilia_Romagna 944689 2109482 795045 2228497 4337979
7 Friuli_Venezia_Giulia 261186 596265 212969 634671 1230936
8 Lazio 1241962 2703994 1126148 2922716 5626710
9 Liguria 320014 767057 273659 848007 1615064
10 Lombardia 2148215 4762370 1806394 4980306 9742676
11 Marche 334044 763741 278521 805837 1569578
12 Molise 68308 156036 58492 164759 320795
13 Piemonte 917932 2149373 761839 2283198 4432571
14 Sardegna 399198 819518 349360 851483 1671001
15 Sicilia 1113415 2433605 1009358 2604194 5037799
16 Toscana 757715 1787668 649483 1920150 3707818
17 Trentino_Alto_Adige 251109 500811 215719 517846 1018657
18 Umbria 182907 431313 155483 462909 894222
19 Valle_D_Aosta 29314 62451 23671 64614 127065
20 Veneto 1083455 2392663 905460 2492885 4885548
Dal sito Istat abbiamo importato questa tabella con la colonna di donne nubili e la colonna di donne totali, per regione. A questo punto calcoliamo la percentuale relativa per ogni regione di donne nubili:
perc_nubili <- data$Femmine_Nubili / data$Totale_Femmine * 100
Dividiamo adesso questo vettore in un fattore a 5 livelli in base alla percentuale (così da ottenere ad esempio le classi: "altissima percentuale", "alta percentuale", "media percentuale", "bassa percentuale", "bassissima percentuale", di donne nubili). Utilizziamo la funzione
cut2 della library(Hmisc):
library(Hmisc)
lev <- cut2(perc_nubili, g=5)
Il parametro
g=5 indica proprio il numero di classi. Automaticamente la funzione calcola i cut-off in modo tale da dividere i 20 valori in 5 gruppi, ciascuno con 4 valori.Ora dobbiamo creare un vettore contenente l'elenco delle classi, relativi alle regioni. Ossia nell'ordine delle regioni visto prima, dobbiamo ottenere un vettore che assegni il livello a ciascuna regione (ad esempio la prima regione, Abruzzo, appartiene alla quarta classe). Utilizziamo la seguente stringa:
nubili <- as.character(as.numeric(lev))
nubili
[1] "2" "4" "4" "4" "5" "3" "1" "4" "1" "3" "2" "2" "1" "5" "5" "2" "5" "1" "3" "3"
Adesso creiamo un fattore dei livelli così ottenuti, e lo aggiungiamo alla variabile gadm, contenente i dati della mappa:
gadm$nubili <- as.factor(nubili)
Creiamo infine un vettore contenente i 5 colori (tanti quanti sono i livelli creati):
col <- rev(heat.colors(length(levels(lev))))
heat.colors permette di ottenere una palette di colori sulle tonalità rosse; altre palette già impostate in R sono:-
rainbow-
heat.colors-
terrain.colors-
topo.colors-
cm.colorsA questo punto possiamo colorare la nostra mappa. Utilizziamo la funzione spplot della library sp. La sintassi che ho utilizzato è la seguente:
-
gadm: è il nome della variabile contenente i dati della mappa-
"nubili": è il vettore contenente i livelli per ciascuna regione-
col.regions: è il vettore contenente i colori da utilizzare-
main: è il titolo della finestra
spplot(gadm, "nubili", col.regions=col, main="Classi regionali di donne nubili")
Il risultato ottenuto è il seguente:
Come si può vedere dalla mappa, le regioni a più bassa percentuale di donne nubili (colore giallo chiaro, classe 1) sono: Umbria, Friuli Venezia Giulia, Piemonte e Liguria; le regioni a più alta percentuale (colore rosso, classe 5) invece sono: Sicilia, Sardegna, Campania, Trentino Alto Adige.
Sul sito Gadm, sono anche disponibili mappe più dettagliate, che dividono l'Italia in base alle province, e anche in base alle città, in modo da ottenere choropleth map ancora più dettagliate.
Per par condicio, ecco il codice e la mappa relativi alla percentuale di maschi celibi:
perc_celibi <- data$Maschi_Celibi / data$Totale_Maschi * 100
data2 <- data.frame(data, perc_celibi)
lev <- cut2(perc_celibi, g=5)
con <- url("http://gadm.org/data/rda/ITA_adm1.RData")
print(load(con))
close(con)
celibi <- as.character(as.numeric(lev))
gadm$celibi <- as.factor(celibi)
col <- rev(heat.colors(length(levels(lev))))
spplot(gadm, "celibi", col.regions=col, main="Classi regionali di uomini celibi")
