Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
R code to assure reproducibility CNR web location data repo http://149.139.8.55/data/comuni_ispra_suolo/ . See ISPRA licence for data reuse.
###########################################################################
# Cutting and reprojecting in WGS84 Geographic the ISPRA data of soil leasing for each Administrative unit - Italy following
# OpenstreetMap Extracts bounds https://github.com/osmItalia/estratti-locali-openstreetmap.
# IBIMET CNR for @CNRconsumosuolo twitter
# Vivaioricerca
##################################################################################################
if (!require("RColorBrewer")) {
install.packages("RColorBrewer")
library(RColorBrewer)
}
if (!require("raster")) {
install.packages("raster")
library((raster)
}
if (!require("rgdal")) {
install.packages("rgdal")
library(raster)
}
if (!require("spatial.tools")) {
install.packages("spatial.tools")
library(spatial.tools)
}
if (!require("leaflet")) {
require(devtools)
devtools::install_github("rstudio/leaflet")
library(leaflet)
}
########################################################################################
load("comuni_ita_bounds.rda")
################################################################################################
ISPRA_proj="+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
newproj=CRS("+init=epsg:4326")
########################################################################################
uso_suolo_ita=raster("imp_5m.tif")
# dimensions : 272103, 214715, 58424595645 (nrow, ncol, ncell)
########################################################################################
for (i in 1: nrow(comuni_ita_bounds)) {
e <- extent(comuni_ita_bounds$minlon[i],
comuni_ita_bounds$maxlon[i],
comuni_ita_bounds$minlat[i],
comuni_ita_bounds$maxlat[i])
box_consumo<-bbox_to_SpatialPolygons(e, proj4string = CRS("+init=epsg:4326"))
box_consumo_ISPRA_proj<- projectExtent(box_consumo, ISPRA_proj)
uso_suolo_ita_cut=crop(uso_suolo_ita,box_consumo_ISPRA_proj)
writeRaster(uso_suolo_ita_cut, filename=paste0(comuni_ita_bounds$name[i],"_ispra_csuolo_or.tif"), format="GTiff", overwrite=TRUE)
system(paste0("gdalwarp -t_srs EPSG:4326 ",comuni_ita_bounds$name[i],"_ispra_csuolo_or.tif ", comuni_ita_bounds$name[i],"_ispra_csuolo.tif"))
unlink(paste0(comuni_ita_bounds$name[i],"_ispra_csuolo_or.tif"))
unlink(paste0(comuni_ita_bounds$name[i],"_ispra_csuolo_or.xml"))
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.