User Tools

Site Tools


wiki:map_rf_itaveg_ens2028_r

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

wiki:map_rf_itaveg_ens2028_r [2017/12/05 22:53] (current)
Line 1: Line 1:
 +====== Mapping the Future projections of distribution for Natural Forest in Italy ======
 +   
 +The objective is to use the randomForest distribution model and a series of 23 environmental predictor variables including future climatic conditions, to map natural forests. \\   
 +The bioclimatic predictors are issued from the IPCC SRES A2a Scenario. The climatic model is an averaged computation of CISRO, CCCMA and HADCM3 models. ​  (see www.worldclim.org for specifications and data download) .\\
 +
 + see /​../​map_RF_ITAveg.R ​ for specifications on the functioning of the present script
 +
 +==== Load required libraries and import randomforest model ====
 +<code r>
 +   ​library(rgdal)
 +   ​library(maptools)
 +   ​library(grid)
 +   ​library(randomForest) ​
 +</​code>​
 +
 +If necessary you will have to install some packages as follows:
 +<code r>
 +install.packages(pkgs="​maptools",​lib= "/​usr/​lib/​R/​site-library/"​)
 +</​code>​
 +
 +Import model Random forest model:
 +<code r>
 +load("​~/​ost4sem/​studycase/​ITA_veg/​models/​RF23tune_ITAveg"​)
 +gc()
 +</​code> ​
 +
 +==== Define your working path ====
 +
 +<code r>
 +  path.rasters = as.character("/​home/​ste/​ost4sem/​grassdbnew/​europe/​Vmodel/​cellhd/"​)
 +  path.results = as.character("/​home/​ste/​ost4sem/​studycase/​ITA_veg/​outmaps/"​)
 +</​code>​
 +
 +
 +
 +==== Create a predictor grid stack ====
 +<code r>
 +   ​print("​ ------------- NOW CREATE predictors grid stack  -------------"​)
 +</​code> ​
 +
 +Importing ARC INFO grid into R SpatialGridDataFrame and generate a spatial data stack with readGDAL function ​
 +This is done in a loop tiles to tiles
 +<code r>
 +   tile = c(20*0:114)
 +   
 +   for(t in 1:57){
 +   
 +   ​datagrids = readGDAL(paste(path.rasters,"​pr104ENSA2a80",​sep=""​),​ offset=c(0,​tile[t]),​ region.dim=c(1480,​20))
 +   ​datagrids$pr104 = datagrids$band1
 +   ​datagrids$pr105 = readGDAL(paste(path.rasters,"​pr105ENSA2a80",​sep=""​),​ offset=c(0,​tile[t]),​ region.dim=c(1480,​20))$band1
 +   
 +   ​datagrids$pr106 = readGDAL(paste(path.rasters,"​pr106ENSA2a80",​sep=""​),​ offset=c(0,​tile[t]),​ region.dim=c(1480,​20))$band1
 +   gc()
 +   
 +   ​datagrids$pr109 = readGDAL(paste(path.rasters,"​pr109ENSA2a80",​sep=""​),​ offset=c(0,​tile[t]),​ region.dim=c(1480,​20))$band1
 +   ​datagrids$pr110 = readGDAL(paste(path.rasters,"​pr110ENSA2a80",​sep=""​),​ offset=c(0,​tile[t]),​ region.dim=c(1480,​20))$band1
 +   
 +   gc()
 +   ​datagrids$pr111 = readGDAL(paste(path.rasters,"​pr111ENSA2a80",​sep=""​),​ offset=c(0,​tile[t]),​ region.dim=c(1480,​20))$band1
 +   ​datagrids$pr118 = readGDAL(paste(path.rasters,"​pr118ENSA2a80",​sep=""​),​ offset=c(0,​tile[t]),​ region.dim=c(1480,​20))$band1
 +   
 +   gc()
 +   ​datagrids$pr120 = readGDAL(paste(path.rasters,"​pr120ENSA2a80",​sep=""​),​ offset=c(0,​tile[t]),​ region.dim=c(1480,​20))$band1
 +   ​datagrids$pr123 = readGDAL(paste(path.rasters,"​pr123ENSA2a80",​sep=""​),​ offset=c(0,​tile[t]),​ region.dim=c(1480,​20))$band1
 +      ​
 +   gc()
 +   ​datagrids$pr126 = readGDAL(paste(path.rasters,"​pr126ENSA2a80",​sep=""​),​ offset=c(0,​tile[t]),​ region.dim=c(1480,​20))$band1
 +   ​datagrids$pr128 = readGDAL(paste(path.rasters,"​pr128ENSA2a80",​sep=""​),​ offset=c(0,​tile[t]),​ region.dim=c(1480,​20))$band1
 +   
 +   gc()
 +   ​datagrids$pr131 = readGDAL(paste(path.rasters,"​pr131ENSA2a80",​sep=""​),​ offset=c(0,​tile[t]),​ region.dim=c(1480,​20))$band1
 +   ​datagrids$pr134 = readGDAL(paste(path.rasters,"​pr134ENSA2a80",​sep=""​),​ offset=c(0,​tile[t]),​ region.dim=c(1480,​20))$band1
 +   
 +   gc()
 +   ​datagrids$pr135 = readGDAL(paste(path.rasters,"​pr135ENSA2a80",​sep=""​),​ offset=c(0,​tile[t]),​ region.dim=c(1480,​20))$band1
 +   ​datagrids$pr136 = readGDAL(paste(path.rasters,"​pr136ENSA2a80",​sep=""​),​ offset=c(0,​tile[t]),​ region.dim=c(1480,​20))$band1
 +   
 +   gc()
 +   ​datagrids$pr137 = readGDAL(paste(path.rasters,"​pr137ENSA2a80",​sep=""​),​ offset=c(0,​tile[t]),​ region.dim=c(1480,​20))$band1
 +   ​datagrids$pr138 = readGDAL(paste(path.rasters,"​pr138ENSA2a80",​sep=""​),​ offset=c(0,​tile[t]),​ region.dim=c(1480,​20))$band1
 +   
 +   ​datagrids$pr198 = readGDAL(paste(path.rasters,"​pr198ENSA2a80",​sep=""​),​ offset=c(0,​tile[t]),​ region.dim=c(1480,​20))$band1
 +   ​datagrids$pr200 = readGDAL(paste(path.rasters,"​pr200ita",​sep=""​),​ offset=c(0,​tile[t]),​ region.dim=c(1480,​20))$band1
 +   gc()
 +   
 +   ​datagrids$pr201 = readGDAL(paste(path.rasters,"​pr201ita",​sep=""​),​ offset=c(0,​tile[t]),​ region.dim=c(1480,​20))$band1
 +   ​datagrids$pr202 = readGDAL(paste(path.rasters,"​pr202ita",​sep=""​),​ offset=c(0,​tile[t]),​ region.dim=c(1480,​20))$band1
 +   gc()
 +   ​datagrids$pr203 = readGDAL(paste(path.rasters,"​pr203ita",​sep=""​),​ offset=c(0,​tile[t]),​ region.dim=c(1480,​20))$band1
 +   ​datagrids$pr204 = readGDAL(paste(path.rasters,"​pr204ita",​sep=""​),​ offset=c(0,​tile[t]),​ region.dim=c(1480,​20))$band1
 +   
 +   gc()
 +   
 +   ​cat("​ ---- 23 grids inported into datagrids file - now drawing the map------\n"​)
 +</​code>​
 +
 +==== Imput future climate into the Random Forest model ====
 +The predict function is used to calculate pixel values and index according to the random forest model and original rasters imported ​
 +<code r>
 +   ​pred.mapRF = predict (RF23tuned , newdata=datagrids,​ type="​response"​)
 +   ​cat("​class predicted\n"​)
 +   
 +   ​pred.mapRFindex = predict (RF23tuned , newdata=datagrids,​ type="​vote"​)
 +   ​cat("​class index recorded\n"​)
 +   
 +   ​cat("​map values for Forest Type calculated and added to predict map vector\n"​)
 +</​code>​
 +
 +Generate an empty tile with geographical coordinates to pixel values =-9999
 +<code r>
 +   ​data2map = readGDAL(paste(path.rasters,"​pr104ita",​sep=""​),​ offset=c(0,​tile[t]),​ region.dim=c(1500,​20))
 +   ​cat("​geocoordinates inported\n"​)
 +   
 +   ​data2map$band1 = 20
 +   ​cat("​ no data value is now = 15 in this tile[t]\n"​)
 +</​code>​
 +
 +Add predicted values to the empty tile[t] ​
 +<code r>
 +   ​data2map$band1[as.numeric(dimnames(pred.mapRFindex)[[1]])] = as.numeric(pred.mapRF)
 +   ​cat("​class values added to tile[t]\n"​)
 +</​code>​
 +
 +Generate ascii map tiles ====
 +Generate an asciii TILE file for modelled species ​
 +<code r>
 +   ​write.asciigrid(data2map,​ paste(path.results,​t,"​tile_ITAveg_ENS80.asc",​sep=""​))
 +   
 +   ​cat("​-----------------------------------------Tile",​ t,"​created--------------------------\n"​)
 +</​code>​
 +
 +
 +==== Clean up the memory and close the tile loop ====
 +<code r>
 +rm(data2map)
 +rm(datagrids)
 +rm(pred.mapRF)
 +rm(pred.mapRFindex)
 +gc()
 +
 +}
 +</​code>​
 +
 +
 +
 +
 +==== Merging tiles ====
 +Now we use a bash shell command from the gdal library to merge tiles into the final mosaic map of modeled current natuaral forest vegetation of Italy and we delete these once the map is produced.
 +   
 +<code r>
 +  system("​cd ost4sem/​studycase/​ITA_veg/​results/​maps"​)
 +  system("​gdalwarp ​ -srcnodata "​-9999" ​ -dstnodata "​-9999"​ -of GTiff -ot Int16 -wt Int16     ​*tile_ITAveg_ENS80.asc ​  ​ITAveg_A2a_ENS80.tif"​) ​
 +  system("​rm *tile_ITAveg.asc"​)
 +</​code>​
 +
 +This script is finished.
 +You have now mapped current and future natural forest vegetation in Italy.
  
wiki/map_rf_itaveg_ens2028_r.txt ยท Last modified: 2017/12/05 22:53 (external edit)