############################################################################## # # MODULE: pfailroc.r # AUTHOR: Martin Mergili # # PURPOSE: Slope stability calculations # with the infinite slope stability model # Random variation of parameters # to compute slope failure probability # Script for the creation of an ROC plot # # COPYRIGHT: (c) 2015 - 2017 by the author # (c) 2013 - 2017 by the BOKU University, Vienna # (c) 2015 - 2017 by the University of Vienna # (c) 1993 - 2017 by the R Development Core Team # # VERSION: 20170210 (10 February 2017) # # This program is free software under the GNU General Public # License. # ############################################################################## #Loading libraries library(sp) library(rgdal) library(gplots) library(ROCR) #Importing arguments args<-commandArgs(trailingOnly = TRUE) outdir<-as.name(args[1]) rpfail<-as.character(args[2]) rinventory<-as.character(args[3]) #Creating plot file rocplot<-paste(outdir, 'roc.png', sep='') png(filename=rocplot, width=8, height=7, units='cm', res=300) #Defining margins par(mar=c(2.5,2.5,0.5,0.5)) #Initializing plot and drawing random prediction line plot(x=c(-2,2), y=c(-2,2), col='gray', lty=3, lwd=1, xlim=c(-0.05,1.0), ylim=c(0.0,1.02), axes=FALSE, xlab=NA, ylab=NA) clip(0,1,0,1) abline(a=0, b=1, col='gray', lty=2, lwd=1) clip(-1,2,-1,2) #Defining graphical parameters lwidth='1' lltp=2 lcol='chartreuse4' ia=paste('', sep='') #Reading input data observed<-readGDAL(fname=rinventory) probability<-readGDAL(fname=rpfail) #Preprocessing data frames probability$band1[probability$band1<0]<-NaN observed$band1[observed$band1>0]<-1 observed$band1[observed$band1<0]<-NaN probability$band1[is.na(observed$band1)==TRUE]<-NaN observed$band1[is.na(probability$band1)==TRUE]<-NaN #Generating ROC curve roc<-prediction(probability$band1, observed$band1) #ROC prediction proc<-performance(roc, 'tpr', 'fpr') #performance aucrocv<-performance(roc, 'auc') #area under curve #Formatting area under curve value for display aucroc<-vector('expression', 1) aucroch<-vector('expression', 1) aucroc[1]<-substitute(expression(iaa~aucroca), list(iaa = format(ia), aucroca = format(round(as.numeric(aucrocv@y.values), digits=3), nsmall=3)))[2] #Plotting ROC curve plot(proc, add=TRUE, axes=F, xlab=NA, ylab=NA, lwd=lwidth, col=lcol, lty=lltp) #Plotting bounding box and axes box() axis(side=1, tck=-0.02, labels=NA, at=c(0.0,0.2,0.4,0.6,0.8,1.0)) axis(side=2, tck=-0.02, labels=NA, at=c(0.0,0.2,0.4,0.6,0.8,1.0)) axis(side=1, lwd=0, line=-0.7, at=c(0.0,0.2,0.4,0.6,0.8,1.0)) axis(side=2, lwd=0, line=-0.7, at=c(0.0,0.2,0.4,0.6,0.8,1.0)) #Plotting axis labels mtext(side=1, expression(italic(r)[FP]/italic(r)[ON]), line=1.3) mtext(side=2, expression(italic(r)[TP]/italic(r)[OP]), line=1.3) #Formatting string for display aucroch<-expression(italic(AUC)[ROC]) #Displaying legend legend('bottomright', pch=NA, lty=c(0,lltp), lwd=c(0,lwidth,lwidth,lwidth,lwidth), col=c('black',lcol), legend=c(aucroch,aucroc), y.intersp=1, text.col=c('black',lcol), bty='n', adj=c(0.1,0.5)) #Closing plot file dev.off()