Un articolo veloce come lo script che allego. Se vi siete imbattuti nel Geoportale della Regione Calabria avrete notato sicuramente la sezione OpenData. Tutti i dati sono scaricabili liberamente e si trovano all’interno di un server FTP. Non è complicato scaricare tutto il necessario utilizzando un client FTP, ma mi piaceva condividere un piccolo script scritto in R che consente di scaricare ed unire (mosaic) in un colpo solo tutti gli elementi del DTM a livello comunale.
Come funziona? Lo script crea una funzione nella sessione di lavoro di R, alla funzione vanno specificati due argomenti: il nome del Comune ed il percorso al quadro di unione. Una volta eseguita viene creata una cartella ‘download’ nella directory corrente all’interno della quale saranno scaricati i file *.img relativi agli elementi del DTM. Gli elementi vengono selezionati sulla base di un’intersezione spaziali tra il limite amministrativo del Comune scelto ed il quadro di unione (link sotto) degli elementi del DTM. Una volta terminato il processo di scarico, gli elementi vengono “mosaicati” per ottenere un unico file in formato GeoTIFF.
Di seguito un esempio di utilizzo della funzione per scaricare ed ottenere un unico file georeferenziato in formato GeoTiff del comune di Badolato (CZ):
> download.dtm.calabria('BADOLATO', "/Users/slarosa/Downloads/Q_UNIONE_05.shp")
semplice vero?
Ecco lo script:
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| download.dtm.calabria <- function(comune, qu) { | |
| require("RCurl", quietly = TRUE) | |
| require(downloader, quietly = TRUE) | |
| require(rgdal, quietly = TRUE) | |
| require(gdalUtils, quietly = TRUE) | |
| if (missing(comune)) | |
| stop("Error: the 'comune' argument is mandatory.") | |
| if (missing(qu)) | |
| stop("Error: the 'qu' argument is mandatory.") | |
| dir.out <- paste0('download/') | |
| dir.create(dir.out, showWarnings = FALSE) | |
| url <- "ftp://ftpopendata:OPENDATA2013@geoportale.regione.calabria.it/DTM5X5/" | |
| wfs.url <- 'WFS:http://wms.pcn.minambiente.it/ogc?map=/ms_ogc/wfs/LimitiAmministrativi_2011.map&version=1.0.0&request=GetFeature&typeName=UA.UNITAAMMINISTRATIVE.COMUNI' | |
| wfs.filter <- paste0( | |
| '&filter=<Filter><PropertyIsEqualTo><PropertyName>comune</PropertyName><Literal>', | |
| comune, | |
| '</Literal></PropertyIsEqualTo></Filter>' | |
| ) | |
| grid <- readOGR(qu) | |
| proj4string(grid) <- CRS('+init=epsg:32633') | |
| ogr2ogr( | |
| paste0(wfs.url, wfs.filter), | |
| file.path(getwd(), dir.out), | |
| "UA.UNITAAMMINISTRATIVE.COMUNI", | |
| f = "ESRI Shapefile" | |
| ) | |
| shp <- readOGR(paste0(dir.out, "UA.UNITAAMMINISTRATIVE.COMUNI.shp"), stringsAsFactors = FALSE) | |
| shp <- spTransform(shp, CRS('+init=epsg:32633')) | |
| grid.intersection <- gIntersection(grid, shp, byid = TRUE) | |
| grid <- grid[grid.intersection, ] | |
| cnt <- 1 | |
| for (i in seq(dim(grid@data)[1])) { | |
| res <- try(download(paste0(url, | |
| "grid", grid@data$elemento[i], '.img'), | |
| paste0(dir.out, "grid", grid@data$elemento[i], '.img'), | |
| mode = "wb")) | |
| if (class(res) == 'try-error') { | |
| print('ERROR') | |
| } | |
| print(paste0(cnt, ' of ', dim(grid@data)[1])) | |
| cnt <- cnt + 1 | |
| } | |
| raster.files <- list.files( | |
| path = dir.out, | |
| full.names = TRUE, | |
| pattern = "^grid", | |
| all.files = FALSE | |
| ) | |
| mosaic_rasters( | |
| gdalfile = raster.files, | |
| dst_dataset = paste0(dir.out, "DTM5M.tif"), | |
| of = "GTiff", | |
| a_srs = '+init=epsg:32633 +proj=utm +zone=33 +datum=WGS84 +units=m +no_defs +ellps=WGS84 ' | |
| ) | |
| tmp.file <- | |
| list.files(path = dir.out, | |
| pattern = "UA.UNITAAMMINISTRATIVE.COMUNI|.img", | |
| full.names = T) | |
| file.remove(tmp.file) | |
| } |
ed un piccolo video:
ed anche lo shapefile relativo al quadro di unione: link
.salvatorelarosa 