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")
Ciao,
RispondiEliminaho una domanda:
che tipo di test devo usare se ho una sola variabile indipendente dicotomica che è ad es. "mutazione del gene G si" o "mutazione no" e tale variabile mi influenza 6 o 7 altre variabili dicotomiche dipendenti tipo "morto si" e "morto no" o "operato si" e "operato no", "metastasi si" "metastasi no" ecc..? Quale analisi multivariata dovrei usare? Che comandi servono per R?
Ripeto: una dicotomica che influenza altre variabili dicotomiche, che multivariata utilizzo? quali comandi con R? Per caso è una loglineare? O un'analisi discriminante? O che cosa?
Andrea
questo blog è formidabile e tu sinceramente sei un grande, complimenti vivissimi.
RispondiElimina@Anonimo:
RispondiEliminafin troppo gentile! grazie :D
@Andrea:
ti ho risposto per e-mail, magari confermami se ti è arrivata correttamente :)
una domanda: come si può fare per usare solo una porzione della mappa gadm in versione R??
RispondiEliminaMi spiego meglio: se ad esempio si desiderano colorare i territori dei comuni veneti (level 3 ) bisogna procurarsi il file relativo (che su gadm non c'è!) oppure c'è un modo semplice per ovviare al problema?
Ciao,
RispondiEliminano, non credo sia possibile quello che hai chiesto, con le mappe gadm. Le uniche soluzioni che per ora mi vengono in mente, sono tuttavia:
1) utilizzare una mappa in formato SVG del territorio veneto, da importare in R (esiste un pacchetto apposito), da colorare con mezzi del tutto analoghi (le mappe SVG per intenderci, sono quelle utilizzate da WikiPedia);
2) utilizzare la mappa di livello 3 di gadm, colorando solamente i territori di tuo interesse. Quindi zoommare la regione veneta, con funzioni di zoom già esistenti (prova a cercare su google inglese: "zoom plot R"). In seguito acquisisci l'immagine ed eventualmente la modifichi con un altro programma, per cancellare i territori contigui (non dovrebbe essere troppo difficile).
Fammi sapere se risolvi il problema :)
Ciao Sono Alessia,
RispondiEliminaho trovato molto interessanto questo post.. e vorrei poterlo applicare a livello provinciale e anche per singolo comune.. ma nn riesco a capire l'ordine con ciu vengono caricati i dati geografici dal sito..
infatti le province all'interno delle regioni nn seguono l'ordine alfabetico..
inoltre mi piacrebbe poter "zoommare" su una zona...
Mi piacerebbe tanto avere un aiuto...
:) Alessia M.
Ciao Alessia,
RispondiEliminaallora, per quanto riguarda l'ordine con cui vengono presentate le province e i comuni, fai così:
a) se ti interessato le province, carica la mappa level2 dal sito GADM (con il codice che ho scritto nel post, modificando il link), e poi digita gadm$NAME_2 - in questo modo otterrai l'elenco delle province, nell'ordine con cui dovrai presentare i tuoi dati.
b) se ti interessano i comuni, carica la mappa level3, e poi digita: gadm$NAME_3
(l'ordine in entrambi in casi è questo: ordine alfabetico delle regione, con sottogruppo in ordine alfabetico delle province, e sottogruppo in ordine alfabetico dei comuni; quindi non c'è un ordine alfabetico totale, ma ramificato).
Per quanto riguarda la zoom, come ho scritto appena sopra, prova a cercare su google inglese "zoom plot R", e troverai alcune mailing-list dove parlano di questo. Vedrai che sarà facile specificare facilmente la porzione del grafico ottenuto da zoommare con alcune funzioni già pre-impostate.
:)
Andrea se hai un database sui tumori che non si trova su internet ti chiedo gentilmente di inviarlo a puccigian@freemail.it specialmente se hai anche(oltre i campi che hai indicato) la data operazione,chemio si e chemio no,data decesso etc.
RispondiEliminaNon chiedo nome e cognome ma solo codici per individuare l'individuo, è una cosa che mi riguarda personalmente.
Ciao,
RispondiEliminaun modo alternativo per creare le mappe è l'utilizzo dei shapefile (.shp).
Per utilizzarli è sufficiente caricarli al posto delle mappe Gadm.
ecco il codice:
al posto di
1. con <- url("http://gadm.org/data/rda/ITA_adm1.RData")
2. print(load(con))
3. close(con)
inserire:
tv_c <- readShapePoly("/Users/harloc/Desktop/geo_R/TV_comuni/tv_comuni.shp")
e seguire passo passo la procedura descritta precedentemente.
Ecco dove trovare i file shp:
http://wiki.gfoss.it/index.php/Geodati_Regioni
Complimenti per il sito e in particolare per questo articolo, che mi ha permesso di approffondire R.
Ciao Marco
Ciao,
RispondiEliminati volevo chiedere aiuto su come posso fare la cartina e differenziare il colore associato a ciascuna regione
ho seguito le tue indicazioni ma mi ha errore quando eseguo il comando
print(load(con))
mi sai dare una spiegazione?
grazie
Francesca
Ciao Francesca,
RispondiEliminahai provato con gli shapefile?
ciao marco
Ciao,
RispondiEliminainnanzitutto volevo farti i miei complimenti per quello che sai fare!sei grande!
Poi volevo chiederti se è possibile aggiungere i nomi a ciascuna regione nel grafico.
Grazie mille in anticipo.
Paolo
complimenti per tutto il tuo meraviglioso blog. Sei geniale e semplice insieme.
RispondiEliminaNon riesco a caricare la mappa bianca dell'Italia, cioè non si apre la formula
con <- url("http://gadm.org/data/rda/ITA_adm1.RData")
print(load(con))
close(con).
Che fare?
Grazie per un eventuale graditissimo aiuto.
Daniele
Ciao Daniele,
RispondiEliminaho appena provato ad usare il solito codice, e riesce a caricare normalmente. Forse è un problema di connessione (riprova tra qualche minuto), oppure qualche firewall blocca il programma sul tuo pc.
Se non dovesse risolversi il problema, incolla qui l'errore che ti viene dato, e vedo se posso darti un aiuto :)
ps grazie per i complimenti!! :)
Ti ringrazio per la risposta rapida. Ho riprovato e purtroppo non riesco. Non riesco a copiare close(con) perchè appena clicco dopo la seconda riga la risposta è
RispondiElimina[1]gadm. Ti copio tutto
il pacchetto 'Hmisc' è stato creato con R versione 2.11.0
> local({pkg <- select.list(sort(.packages(all.available = TRUE)))
+ if(nchar(pkg)) library(pkg, character.only=TRUE)})
Warning message:
il pacchetto 'sp' è stato creato con R versione 2.11.0
> # con <- url("http://gadm.org/data/rda/ITA_adm1.RData")
> con <- url("http://gadm.org/data/rda/ITA_adm1.RData")
> print(load(con))
[1] "gadm"
> close(con)
Ancora grazie
Daniele
E' corretto, deve uscire [1] "gadm", dopo il comando print(). Non c'è nell'articolo, perchè in genere quando scrivo il codice, non lo copio dal programma, altrimenti sarebbe comparso.
RispondiEliminaPer quel che vedo non ci sono errori (apparte l'aggiornamento del programma che ti consiglio di fare), quindi credo che puoi procedere tranquillamente senza incontrare altri problemi.
Tu sei veramente un grande. Non ho parole per esprimere la gratitudine per il tuo lavoro. Sono riuscito nell'intento di disegnare la mappa bianca dell'italia digitando Plot(gadm) e miracolosamente è apparsa la cartina bianca dell'Italia.
RispondiEliminaIl tuo lavoro deve essere conosciuto e apprezzato, non sono molti quelli che possiedono la tua chiarezza, sicuramente nemmeno in ambito universitario.
Grazie ancora per il tuo tempo
Daniele (daxgi09@gmail.com)
Ciao,
RispondiEliminain riferimento alla tua risposta del 13 gennaio 2010, ho provato tutti e due i tentativi che tu proponi ma nessuno dei due sembra funzionare. Per essere precisi, io sto cercando di plottare una sola regione italiana con tutti i comuni e di colorare la mappa in base a dei livelli dati.
Ti spiego di seguito i problemi che ho incontrato.
1)Su wikipedia non sono presenti le mappe delle sole regioni con tutti i comuni ma ci sono solo le regioni con le province;
2)la mappe in gadm per essere visualizzate necessitano della funzione "spplot" mentre la funzione di zoom sembra funzionare solo con "plot" e quindi non sembra funzionare con le mappe di gadm per nessun livello.
Riusciresti a suggerirmi una funzione per zommare su questo tipo di grafici?Qualcuno è riuscito a risolvere questo problema ossia visualizzare solo una regione con tutti i suoi comuni utilizzando mappe gadm?
Grazie mille a te che hai un sito spettacolare e grazie a tutti!
ciao a tutti,volevo chiedere un consiglio,ho fatto un copia incolla all'interno del programma R del codice scritto qui sopra però non riesco a visualizzare il grafico,come devo fare?devo caricare qualcosa?grazie
RispondiEliminaCiao mitico.
RispondiEliminaSe nella legenda voglio inserire le classi (es. 10-20; 20-30 ecc) invece che i livelli (es. 1; 2 ecc) come faccio?
Ciao e grazie delle preziose info che condividi con tutti.
Ciao, scusa ... ma sono proprio alle prime armi ... e dovrei realizzare una cartina dell'Italia con alcune città evidenziate da bubble di diversi colore/diametri. Mi potresti per favore aiutare?
RispondiElimina1. ho installato R sul mio pc
2. ho scelto il CRAN mirror (Italy Milan)
3. ho installato il pacchetto: maptools
4. ho caricato il pacchetto maptools
5. mi sono creata un file dati.csv con le coordinate interessate ad essere evidenziate
poi ... non sò cosa fare e le ricerche in rete non mi aiutano a risolvere.
mi daresti gentilmente un'indicazione su come procedere?
Ciao,
Annamaria
Ciao, io vorrei inserire delle classi scelte da me invece dei quantili della variabile, ma ho un problema per modificare i colori della cartina:
RispondiElimina> lev <- cut2(dataset$a, cuts=c(8000,13500,19000,24500,30000))
> ms <- as.character(as.numeric(lev))
> gadm$ms <- as.factor(ms)
> col <- rev(terrain.colors(length(levels(lev))))
> spplot(gadm, "ms", col.regions=col)
in quest'ultimo comando mi compare il messaggio di errore:
"Errore in .local(obj, ...) :
length of col.regions should match number of factor levels"
ma il numero dei colori delle regioni è 5, uguale al numero delle classi..
Grazie in anticipo, ciao
ciao.
RispondiEliminaho fatto una cosa simile qui:
http://taffey6977.blogspot.it/2012/11/mappe-comunali-con-r.html
ciao, sto facendo un esame di statistica all uni e ho bisogno,di una mano. Ho due variabili Fertilità e Reddito pro capite delle Nazioni unite e ho bisogno di fare questo grafico per rappresentare la situazione con R, ma il comando spplot non lo legge puoi aiutarmi????
RispondiEliminaCiao, grazie mille per il post. Purtroppo non riesco a caricare la mappa su R dal sito gadm. L'errore è pagina non trovata. Puoi provare anche tu ad eseguire il comando che hai indicato su questo post e verificare che la mappa carica in modo corretto?
RispondiElimina