Esempi di visualizzazioni statistiche su dati GIS con R

R è un software di statistica open-source, con un numero di metodologie di analisi e visualizzazione ben superiore a quello tipico dei software GIS.

Un apparente svantaggio di R è l'uso di righe di comando piuttosto che di una interfaccia grafica. Ma questo consente di copiare e riutilizzare più volte in R una successione di comandi anche complessa.

Anche se è possibile usufruire di R all'interno di Python, Grass e altri, la maniera più semplice, soprattutto per chi si avvicina la prima volta ad R, è quella di utilizzare di R in maniera a se stante.

Un documento di buon livello e free su R è An Introduction to R: Software for Statistical Modelling & Computing di Kuhnert e Venables. Su temi di analisi spaziale, un testo avanzato è Applied Spatial Analysis with R di Bivand, Pebesma e Gómez-Rubio. Per un quadro generale  dei packages esistenti per l'analisi spaziale in R è utile  il sunto che si trova in CRAN Task View: Analysis of Spatial Data.

Per la lettura di dati vettoriali possono essere usati i packages maptools, shapefiles, foreign, RArcInfo, mentre per i raster i pacchetti Rgdal e raster. Esistono anche dei packages specifici per Saga e Grass. Un pacchetto per la geostatistica è gstat.

Lo scopo di questo post è di presentare come sia possibile, in maniera semplice, analizzare dati GIS. Usiamo quindi i packages che risultano più immediati.

Dati vettoriali

In generale i dati vettoriali sono conservati come shapefile o possono essere facilmente convertiti in tale formato con i GIS. A quel punto è facile importarli in R, o leggendo l'intero shapefile, oppure il solo file dbf dello shapefile, nel quale sono conservati gli attributi descrittivi.

Nel primo caso è possibile anche plottare mappe ma questi utilizzi sono quasi sempre eseguibili nei GIS, che invece mancano di funzionalità statistiche avanzate.

Consideriamo quindi il più semplice package foreign, che consente di leggere un qualsiasi file dbf come quelli degli shapefile.

Useremo lo shapefile dei comuni italiani del 2001 (com2001_s.shp), rilasciato dall'ISTAT, come dataset di esempio. A questo shapefile, contenente l'informazione della popolazione, sono stati aggiunti un campo col nome della regione di appartenenza, oltre a due campi che contengono area e perimetro dei comuni.

Dopo avere installato in R il package foreign, lo si carica con il comando:

library(foreign)

Possiamo quindi leggere il nostro file dbf:

comuni <- read.dbf("com2001.dbf")

avendo se necessario l'accortezza di spostarsi nella cartella che contiene lo shapefile (File → Cambia directory).

A questo punto è utile chiamare una funzione di resoconto sui dati letti, per verificare nomi, tipo e statistiche dei campi importati:

summary(comuni)

Verrà visualizzato un elenco di elenco dei campi, con le loro statistiche.

Quelli che utilizzeremo sono i campi POP2001, AREA e COD_REG.

I codici delle regioni sono nella tabella successiva.

1	PIEMONTE
2	VALLE D'AOSTA
3	LOMBARDIA
4	TRENTINO-ALTO ADIGE
5	VENETO
6	FRIULI-VENEZIA GIULIA
7	LIGURIA
8	EMILIA-ROMAGNA
9	TOSCANA
10	UMBRIA
11	MARCHE
12	LAZIO
13	ABRUZZO
14	MOLISE
15	CAMPANIA
16	PUGLIA
17	BASILICATA
18	CALABRIA
19	SICILIA
20	SARDEGNA

I singoli campi possono essere richiamati inserendo, fra il nome del dataset e quello dello specifico campo (come compare nella lista prodotta dal comando summary), il simbolo '$'.

Per esempio:

summary(comuni$AREA)

A questo punto possiamo creare grafici di vario tipo, oltre che effettuare analisi statistiche. Seguono alcuni esempi di grafici.

Proviamo a sintetizzare con boxplot le distribuzioni dei residenti per comune, raggruppati in base alla regione (espressa dal codice Istat).

boxplot(comuni$POP2001 ~ comuni$COD_REG, main='Comuni italiani - popolazione 2001', xlab='Regioni (codice Istat)', ylab='Residenti')

Il codice (semplificato) POP2001 ~ COD_REG esprime la suddivisione dei valori di popolazione POP2001 in base a valori con significato di categoria conservati in COD_REG.

La statistica è dominata dalle grandi città (Roma, Milano, Napoli, Torino etc.), con le loro popolazioni sull'ordine del milione di abitanti. La grande maggioranza dei comuni finisce schiacciata verso il basso, per il loro numero molto ridotto di residenti.

Per rendersi conto della situazione dei comuni medio-piccoli, è quindi necessario restringere il dataset, con una semplice operazione di filtro sui comuni con meno di 100.000 abitanti ([comuni$POP2001 < 100000]) che va ripetuta per ogni campo coinvolto (due in questo caso). 

boxplot(comuni$POP2001[comuni$POP2001 < 100000] ~ comuni$COD_REG[comuni$POP2001 < 100000], col='pink', main='Comuni italiani con popolazione 2001 < 100.000 residenti', xlab='Regioni (codice Istat)', ylab='Residenti')

Dal boxplot si nota come Puglia (16) e Toscana (9) abbiano distribuzioni meno frammentate rispetto alle altre regioni.

Proviamo a comparare tramite istogrammi la densità di due differenti regioni che presentano una differente distribuzione delle popolazioni dei comuni, il Piemonte (1) e la Puglia (16), considerando l'intervallo inferiore a 200,000 residenti.

I due grafici verrano accolti in una singola una finestra ad una riga e 2 colonne usando il comando par(mfrow=c(1,2)).

pop_piemonte <- comuni$POP2001[comuni$COD_REG == 1 & comuni$POP2001 < 200000]

pop_puglia <- comuni$POP2001[comuni$COD_REG == 16 & comuni$POP2001 < 200000]

par(mfrow=c(1,2))

hist(pop_piemonte, col = 'pink', xlim=c(0,2e5), breaks=10, main='Comuni Piemonte < 200,000 res.', xlab='Residenti', ylab='Numero')

hist(pop_puglia, col = 'pink', xlim=c(0,2e5), breaks=20, main='Comuni Puglia < 200,000 res.', xlab='Residenti', ylab='Numero')

Si nota l'abbondante numero di comuni con residenti inferiori a 10,000 abitanti nel caso del Piemonte oltre alla bassa frequenza di comuni con popolazione maggiore di 20,000 abitanti. La situazione è nettamente più equilibrata in Puglia, oltre a corrispondere ad un numero ben inferiore di comuni, di circa un ordine di grandezza.

Per confronto delle aree degli stessi comuni, possiamo utilizzare i seguenti comandi:

aree_piemonte <- comuni$AREA[comuni$COD_REG == 1 & comuni$POP2001 < 200000]

aree_puglia <- comuni$AREA[comuni$COD_REG == 16 & comuni$POP2001 < 200000]

par(mfrow=c(1,2))

hist(aree_piemonte/1e6, breaks=10, xlim=c(0,6e2), col='orange', main='Comuni Piemonte < 200,000 res.', xlab='Area (kmq)', ylab='Numero')

hist(aree_puglia/1e6, breaks=30, xlim=c(0,6e2), col='orange', main='Comuni Puglia < 200,000 res.', xlab='Area (kmq)', ylab='Numero')

Si conferma anche per le aree comunali (con residenti < 200k) lo squilibrio a favore di una regione come il Piemonte rispetto all'esempio della Puglia: il Piemonte presenta una frammentazione territoriale di un ordine di grandezza superiore rispetto alla situazione della Puglia. Occorre tenere presente che le superfici di Piemonte e Puglia sono circa comparabili (19,600 kmq la Puglia, 25,400 kmq il Piemonte).

Dati raster

La struttura dei dati raster è semplice e possono essere esportati come dati tabulari, importabili in R con un comando come read.dbf o read.csv.

Se invece si preferisce importarli come livelli raster, sono utilizzabili packages come rgdal e raster.

Rgdal potrebbe richiedere l'installazione aggiuntiva di GDAL, se non è già presente perché già installata da altri software GIS.

Negli esempi che seguono useremo i dati SRTM2 a 3 secondi-arco, collegati fra loro e poi ritagliati all'estensione del Piemonte. Questi dati presentano alcune zone con dati mancanti.

library(rgdal)

e quindi si possono leggere i formati raster riconosciuti da gdal, come l' ESRI ASCII grid:

piem_elev <- readGDAL("srtm2_orig_geog.asc")

ed ottenere informazioni di sintesi con:

summary(piem_elev)

o visualizzarlo come mappa:

image(piem_elev, col = terrain.colors(12))

Per leggere ESRI ASCII grid esiste anche la funzione readAsciiGrid del package maptools.

Nel formato ESRI grid è presente una sola banda. Possiamo visualizzare un istogramma della sua distribuzione con:

hist(piem_elev$band1, col='orange')

oppure la sua empirical cumulative distribution function (ECDF) con:

plot(ecdf(piem_elev$band1), main='ECDF elevazioni Piemonte', xlab='elevazione')

Questo è solo un breve esempio delle numerose, possibili visualizzazioni con R.