############################################################################## # # MODULE: slopestability.map.r # # AUTHOR: Martin Mergili # # PURPOSE: R script for the preparation # of maps visualizing slope stability model results # # COPYRIGHT: (c) 2013 - 2016 by the author # (c) 2013 - 2016 by the BOKU University, Vienna # (c) 2015 - 2016 by the University of Vienna # (C) 1993 - 2016 by the R Development Core Team # # VERSION: 20160929 (29 September 2016) # # This program is free software under the GNU General Public # License (>=v2). # ############################################################################## #Loading libraries library(sp) library(rgdal) library(raster) #Importing arguments defined in r.slope.stability.py args<-commandArgs(trailingOnly = TRUE) indir<-as.character(args[1]) outdir<-as.name(args[2]) rlayer<-as.character(args[3]) rhillshade<-as.character(args[4]) vinventory<-as.character(args[5]) pflag<-as.integer(args[6]) #Defining titles of maps, color scheme for layers #and custom threshold levels for legend if ( pflag == 0 ) { titlong<-'Factor of safety' titshort<-' FoS' lev0<-0.5 lev1<-0.75 lev2<-1 lev3<-1.3 lev4<-1.6 lev5<-2 lev6<-2.5 col1<-rgb(0.0,0.0,0.0,1.0) col2<-rgb(0.4,0.0,0.4,0.8) col3<-rgb(0.1,0.1,0.8,0.6) col4<-rgb(0.2,0.5,0.8,0.4) col5<-rgb(1.0,1.0,0.0,0.2) col6<-rgb(1.0,1.0,1.0,0.0) } else { titlong<-'Slope failure probability' titshort<-substitute(expression(~~~~~~~~~~~italic(P)[f])) lev0<-0 lev1<-0.1 lev2<-0.2 lev3<-0.4 lev4<-0.6 lev5<-0.8 lev6<-1 col1<-rgb(0.3,0.7,1.0,0.25) col2<-rgb(0.2,0.7,0.6,0.40) col3<-rgb(0.2,0.6,0.2,0.55) col4<-rgb(0.4,0.4,0.2,0.70) col5<-rgb(0.3,0.1,0.1,0.85) col6<-rgb(0.0,0.0,0.0,1.00) } #Defining precision of legend labels rd<-2 #Defining legend labels lab0<-'0.00' lab1<-format(round(lev1,rd),nsmall=rd) lab2<-format(round(lev2,rd),nsmall=rd) lab3<-format(round(lev3,rd),nsmall=rd) lab4<-format(round(lev4,rd),nsmall=rd) lab5<-format(round(lev5,rd),nsmall=rd) if ( pflag == 0 ) { lab6<-'INF' } else { lab6<-'1.00' } #Defining names of raster maps layermap<-rlayer hillshmap<-rhillshade #Importing raster maps layer<-raster(layermap) hillshade<-raster(hillshmap) #Importing vector maps observedshape<-readOGR(indir, layer=vinventory) #Computing x and y extent of raster layers xrmin<-xmin(layer) xrmax<-xmax(layer) yrmin<-ymin(layer) yrmax<-ymax(layer) xdiff<-xrmax-xrmin ydiff<-yrmax-yrmin #Adapting boundaries for vector layers xsmin<-xrmin+0.035*xdiff xsmax<-xrmax-0.035*xdiff ysmin<-yrmin+0.035*ydiff ysmax<-yrmax-0.035*ydiff #Layouting map asprat<-xdiff/(1.4*ydiff) margx<-6.0 margy<-2.5 dispx<-15 dispy<-(15-margx)/asprat+margy #Adapting vector maps if ( pflag == 0 ) { layer[layer==9999]<-NA layer[layer<0.5]<-0.5 layer[layer>2.5]<-2.5 } else if ( pflag == 1 ) { layer[layer==-9999]<-NA layer[layer==0]<-NA } layer[layer>lev6]<-NA layer[layer xsmin ) xaxlbs[j]<-xlmin+(j-1)*ldiff #tick mark position if ( xlmin+(j-1)*ldiff > xsmax-ldiff ) break #break if last label was reached j<-j+1 } j<-1 #initializing counter for loop over all y axis ticks and labels repeat { #starting loop if ( ylmin+(j-1)*ldiff > ysmin ) yaxlbs[j]<-ylmin+(j-1)*ldiff #tick mark position if ( ylmin+(j-1)*ldiff > ysmax-ldiff ) break #break if last label was reached j<-j+1 } #Plotting bounding box and axes box() par(cex=0.65) axis(side=1, tck=-0.01, labels=NA, at=xaxlbs) axis(side=2, tck=-0.01, labels=NA, at=yaxlbs) axis(side=1, lwd=0, line=-0.6, at=xaxlbs) axis(side=2, lwd=0, line=-0.6, at=yaxlbs) #Plotting legend par(cex=1.0) legend('bottomleft', lty=1, lwd=1, col='red', legend='observed landslide areas', text.col='black', bty='n', horiz=TRUE, x.intersp=0.25, inset=c(0.25,0.01)) #Closing plot file dev.off()