Tutorial N. 3 – R a servizio del GIS: Iterare operazioni di definizione e proiezione
Patrocinio Ordine Dottori Agronomi & Forestali delle Prov. di Campobasso e Isernia
27 Gennaio 2018Tutorial N. 4 – Integrare R con QGIS: un matrimonio quasi perfetto!
21 Febbraio 2018Forse non tutti conoscono il linguaggio di programmazione R dato che molto spesso il suo utilizzo rimane confinato in ambito accademico. R nasce soprattutto per l’analisi statistica dei dati, ma il crescente numero di “pacchetti” sviluppati per una moltitudine di scopi, lo sta rendendo sempre più completo aprendo nuovi orizzonti per il suo utilizzo.
Senza ombra di dubbio una delle integrazioni che ha molte potenzialità è quella nel settore GIS. Una serie di pacchetti permette a R di diventare un vero e proprio GIS, con il vantaggio di integrare le funzionalità tipiche dei GIS con analisi statistiche, operazioni complesse, ecc…
Facciamo un piccolo esempio…
Abbiamo una serie di shapefile (centinaia) suddivisi in sottocartelle ai quali manca il file .prj e che devono essere proiettati nel SR EPSG32633. I file sono nel sistema EPSG3004 ma manca il file di proiezione.
Come fare? A mano un per volta? Con un batch processing? (troppo scomodo e macchinoso)
Ecco che ci viene in aiuto il nostro carissimo R (vi consiglio di installare oltre a R anche l’interfaccia RStudio che rende l’approccio più “semplice”).
Apriamo la nostra interfaccia RStudio e iniziamo con installare i pacchetti necessari: il pacchetto “raster” e il pachetto “rgdal” (il nome dovrebbe suonarvi familiare).
library(rgdal)
Settiamo la working directory ovvero la cartella nella quale abbiamo i nostri shapefile
Nella cartella Carta tecnica regionale 5000 ci sono alcune sottocartelle, all’interno delle quali ci sono gli shapefile
Dobbiamo cercare e “listare” tutti gli shapefile presenti nelle sottocartelle tramite il comando list.files
Adesso dobbiamo inserire un ciclo for che permette di applicare una data porzione di codice a più elementi (tutti gli shapefile listati nell’oggetto files)
{
q<-readOGR(files[i])
crs(q)<-“+init=EPSG:3004”
q<-spTransform(q, CRS(“+init=EPSG:32633”))
writeOGR(q, layer=paste(gsub(“^.*\\/”, “”, (gsub(“.shp”, “”, files[i]))), sep = “”),”C:\\Users\\Ludovico\\Desktop\\Carta tecnica regionale 5000\\Output”,driver=”ESRI Shapefile”)
}
- for (i in 1:length(files)) applica la funzione compresa tra { } un numero di volte che varia da 1 a length(files) ovvero al numero di shapefiles contenuti in files. L’indice i assume di volta in volta valori diversi da 1 a length(files)
- q<-readOGR(files[i]) importa lo shapefile i-esimo e lo assegna ad un oggetto generico q
- crs(q)<-“+init=EPSG:3004” assegna la proeizione 3004 allo shapefile i-esimo
- q<-spTransform(q, CRS(“+init=ESPG:32633”)) applica la funzione spTransform all’oggetto q
- writeOGR salva l’oggetto q come shapefile; il nome assegnato è preso dal nome originale dello shapefile tramite la funzione paste;
- paste copia il nome dello shapfile i-esimo completo di percorso: la funzione writeOGR chiede il nome dello shapefile senza estensione e per fare questo combiniamo paste con gsub che permette di eliminare parti di testo
A questo punto non ci resta che premere su Run e lanciare il codice… Andiamo a preparaci un bel caffè e al nostro ritorno troveremo la cartella Output con tutti gli shapefile proiettati in automatico.
Il codice completo
files<-list.files(pattern=”.shp”, recursive=TRUE)
library(raster)
library(rgdal)
for (i in 1:length(files))
{
q<-readOGR(files[i])
crs(q)<-“+init=EPSG:3004”
q<-spTransform(q, CRS(“+init=EPSG:32633”))
writeOGR(q, layer=paste(gsub(“^.*\\/”, “”, (gsub(“.shp”, “”, files[i]))), sep = “”),”C:\\Users\\Ludovico\\Desktop\\Carta tecnica regionale 5000\\Output”,driver=”ESRI Shapefile”)
}
