#*********************************************************************************** #*********************************************************************************** #**************** ********************** #*********** LANDSLIDE VOLUME FIT ***************** #*********** This script was prepared using R 2.6.2 ***************** #**************** by ******************** #************** MariaTeresa.Brunetti@irpi.cnr.it ****************** #**************** ********************** #*********************************************************************************** #*********************************************************************************** # # The self-explanatory R script performs the analysis of landslide volumes through # the following 6 steps: # 1. Read the raw landslide volume data (VL), in linear coordinates (e.g. cubic meters), # from separate text files. # 2. To evaluate the error associated with the scaling exponent (beta), for each dataset, # (i) construct a set of n synthetic datasets, by substituting the empirical measurements # of VL with synthetic values obtained through random sampling from a uniform distribution # spanning a pre-defined volume range, (ii) determine the mean and the standard deviation # of beta. # 3. Transform the raw landslide volume data in logarithmic (base 10) coordinates, log(VL). # 4. Using the 'density' function in R, perform kernel density estimation (KDE) on the # log-transformed landslide volume data, to obtain the probability density p(logVL). # 5. Using Eq. (5) in the text, transform p(logVL) in p(VL), and show the result in a plot # portraying p(VL) as function of VL. # 6. Using the 'nls' function in R, determine the power law fit of the (right) tail of the # log of probability density p(VL) for landslide volumes larger than a threshold value. # #************** set home directory ************** current.dir<-getwd() setwd(current.dir) #************** load "MASS" package to fit statistical distributions ************* library(MASS) #************** STEP 1. ************* name.series<-c("Umbria-Marche","Balza Tagliata") # !!! Change the following paths to your working directory !!! data.table1<-read.table("data/Data_volume_m3_UmbriaMarche.txt", header = TRUE,dec=".", sep="\t") data.table2<-read.table("data/Data_volume_m3_BalzaTagliata.txt", header = TRUE,dec=".", sep="\t") #************** define some parameters ************* kernel.bandwidth<-c(0.4,0.4) start.fit<-c(-1,-1) stop.fit<-c(2.5,2) volume.percent.error<-c(.5,.5) volume.factor.error<-c(2,2) color.id<-c("dark violet","violet") loop.number<-1 #************** define two matrices ************* lin.volume.values<-matrix(nrow=2,ncol=3) colnames(lin.volume.values)=c("Min Value","Max Value","Median Value") rownames(lin.volume.values)<-name.series extreme.volume.values<-matrix(nrow=2,ncol=2) colnames(extreme.volume.values)=c("Left volume limit","Right volume limit") rownames(extreme.volume.values)<-name.series #************** start loop on events ************* for (loop in 1:loop.number) { #************** start for cycle ****************** for (volume.series in 1:2) { series.name<-as.numeric(get(paste("data.table",volume.series,sep=""))[,1]) no.events<-length(series.name) assign(paste("no.events.",volume.series,sep=""),no.events) #*************** Computes the minimum, maximum and the median of the volume series *********** min.value<-min(series.name) lin.volume.values[volume.series,1]<-min.value max.value<-max(series.name) lin.volume.values[volume.series,2]<-max.value median.value<-median(series.name) lin.volume.values[volume.series,3]<-median.value #*************** Computes the minimum, maximum of the volume series with the associate error *********** min.volume.values<-series.name/volume.factor.error[volume.series] max.volume.values<-series.name*(volume.factor.error[volume.series]-volume.percent.error[volume.series]) extreme.volume.values[volume.series,1]<-min(min.volume.values) extreme.volume.values[volume.series,2]<-max(max.volume.values) #************** STEP 2. ************* series.name.random<-runif(no.events,min.volume.values,max.volume.values) #************** STEP 3. ************* series.name<-log10(series.name.random) assign(paste("data.series.volume.",volume.series,sep=""),series.name) total.no.value<-sum(series.name) assign(paste("total.no.value.",volume.series,sep=""),total.no.value) as.numeric(series.name) #************** STEP 4. ************* den<- density(series.name, bw=kernel.bandwidth[volume.series], adjust=1, kernel=c("g"), n=10000, give.Rkern=FALSE) assign(paste("den.",volume.series,sep=""),den) #*********** cut those density values outside the data interval ************* den.selected.values<-which(den$x>=series.name[1] & den$x<=series.name[no.events]) den$x<-den$x[den.selected.values] den$y<-den$y[den.selected.values] #************** STEP 5. ************* kernel.diff<-10^(den$x[2:length(den$x)])-10^(den$x[1:(length(den$x)-1)]) assign(paste("kernel.diff.",volume.series,sep=""),kernel.diff) den.linear<-den$y[1:(length(den$y)-1)]*((den$x[2:length(den$x)])-(den$x[1:(length(den$x)-1)]))/kernel.diff assign(paste("den.linear.",volume.series,sep=""),den.linear) freq.linear<-den$y[1:(length(den$y)-1)]*no.events*((den$x[2:length(den$x)])-(den$x[1:(length(den$x)-1)]))/kernel.diff assign(paste("freq.linear.",volume.series,sep=""),freq.linear) linear.x<-10^(den$x[1:(length(den$x)-1)]) assign(paste("linear.x.",volume.series,sep=""),linear.x) sum(kernel.diff*den.linear) # check for normalization #*********** transform probability density and frequency density in log to fit ************* den.log<-log10(den.linear) assign(paste("den.log.",volume.series,sep=""),den.log) freq.log<-log10(freq.linear) assign(paste("freq.log.",volume.series,sep=""),freq.log) log.x<-log10(linear.x) assign(paste("log.x.",volume.series,sep=""),log.x) #************** STEP 6. ************* nls.control(minFactor = 1/2048) fit.straight.line <- nls(den.log ~ coefficient + (slope*log.x), start = c("coefficient" = 1, "slope" = -1)) fit.straight.line.freq <- nls(freq.log ~ coefficient + (slope*log.x), start = c("coefficient" = 1, "slope" = -1)) value.straight.line.den.theo<-predict(fit.straight.line) value.straight.line.freq.theo<-predict(fit.straight.line.freq) assign(paste("value.straight.line.den.theo.",volume.series,sep=""),value.straight.line.den.theo) assign(paste("value.straight.line.freq.theo.",volume.series,sep=""),value.straight.line.freq.theo) assign(paste("coef.fit.straight.line.",volume.series,sep=""),coef(fit.straight.line)) assign(paste("coef.fit.straight.line.freq.",volume.series,sep=""),coef(fit.straight.line.freq)) #*********** select the interval for the fit ************* data.subset<-as.data.frame(cbind(log.x,freq.log,den.log)) log.x.high<-log.x[log.x>=start.fit[volume.series] & log.x<=stop.fit[volume.series]] data.subset.den.high<-as.data.frame(subset(data.subset, log.x>=start.fit[volume.series]& log.x<=stop.fit[volume.series], select=c(log.x,den.log))) data.subset.freq.high<-as.data.frame(subset(data.subset, log.x>=start.fit[volume.series]& log.x<=stop.fit[volume.series], select=c(log.x,freq.log))) nls.control(minFactor = 1/2048) fit.straight.line.high <- nls(den.log ~ coefficient + (slope*log.x), data.subset.den.high, start = c("coefficient" = -1, "slope" = -1)) fit.straight.line.freq.high <- nls(freq.log ~ coefficient + (slope*log.x), data.subset.freq.high, start = c("coefficient" = 1, "slope" = -1)) value.straight.line.den.theo.high<-predict(fit.straight.line.high) value.straight.line.freq.theo.high<-predict(fit.straight.line.freq.high) assign(paste("log.x.high.",volume.series,sep=""),log.x.high) assign(paste("value.straight.line.den.theo.high.",volume.series,sep=""),value.straight.line.den.theo.high) assign(paste("value.straight.line.freq.theo.high.",volume.series,sep=""),value.straight.line.freq.theo.high) assign(paste("coef.fit.straight.line.high.",volume.series,sep=""),coef(fit.straight.line.high)) assign(paste("coef.fit.straight.line.freq.high.",volume.series,sep=""),coef(fit.straight.line.freq.high)) } #************** end for cycle ****************** pdf(file = paste("Output", loop, ".pdf", sep=""), width = 6, height = 6, onefile = TRUE, family = "Helvetica", fonts = NULL, paper = "special", pagecentre=TRUE) #*********** plot the kernel density of log(VL) ************* # windows() par(mfrow=c(2,1), pty="s", font=6, font.main=6, cex.lab=1.2, cex.main=1.4) for (volume.series in 1:2) { plot(get(paste("den.", volume.series,sep="")), type='l', main=expression(paste("Kernel density of log (", V[L],") - Fig. 1", sep="")), xlab=expression(paste("Log (",V[L],")",sep="")), ylab="Kernel density") mtext(paste(name.series[volume.series]), side=3, col="red", cex=1, line=0.5) rug(get(paste("data.series.volume.", volume.series, sep=""))) } #*********** plot the comparison of kernel densities of log(VL) ************* # windows() par(mfrow=c(1,1), pty="s", font=6, font.main=6, cex.lab=1.2, cex.main=1.4) plot(den.1, type="l", main=expression(paste("Comparison of kernel densities of log (",V[L],") - Fig. 2", sep="")), pch=21, lwd=2, cex=1, col=color.id[1], xlim=c(-5,4), ylim=c(0,0.5), xlab=expression(paste("Log (",V[L],")",sep="")), ylab="Kernel density") legend(0,.5, name.series[1], lwd=2, bty="n", col=color.id[1]) grid(nx = NULL, ny = NULL, col = "lightgray", lty = "dotted", lwd = par("lwd")) lines(get(paste("den.",2,sep="")), type="l", lty=1, pch=21, lwd=2, cex=1, col=color.id[2]) legend(0,.45, name.series[2], lwd=2, lty=1, bty="n", col=color.id[2]) #*********** plot the comparison of frequency densities of log(VL) ************* # windows() par(pty="s", font=6, font.main=6, cex.lab=1.2, cex.main=1.3) plot(den.1$x,den.1$y*no.events.1,type="l",main=expression(paste("Comparison of frequency densities of log", (V[L])," - Fig. 3", sep="")), pch=21, lwd=2, cex=1, col=color.id[1], xlim=c(-5,4), ylim=c(0,50), xlab=expression(paste("Log (", V[L],")",sep="")), ylab="Frequency density") legend(0,250, name.series[1], lwd=2, bty="n", col=color.id[1]) grid(nx = NULL, ny = NULL, col = "lightgray", lty = "dotted", lwd = par("lwd")) lines(den.2$x,den.2$y*no.events.2,type="l", lty=1, pch=21,lwd=2,cex=1,col=color.id[2]) legend(0,235, name.series[2], lwd=2, bty="n", col=color.id[2]) #********* plot the probability density of VL ************* # windows() par(pty="s", font=6, font.main=6, cex.lab=1.2, cex.main=1.4) plot(linear.x.1, den.linear.1, log="xy", type="l", pch=21, lwd=2, cex=1, col=color.id[1], main="Comparison of probability densities - Fig. 4", ylab=expression(paste("p, ( ",m^{-3},")",sep="")), xlab=expression(paste(V[L], ", ( ", m^{3},")",sep="")), xlim=c(10^-5,10^3), ylim=c(10^-6,10^3)) legend(.1,10^2.5, name.series[1], lwd=2, bty="n", col=color.id[1]) grid(nx = NULL, ny = NULL, col = "lightgray", lty = "dotted", lwd = par("lwd")) lines(linear.x.2, den.linear.2, type="l", pch=21, lwd=2, cex=1, col=color.id[2]) legend(.1,10^2, name.series[2], lwd=2, bty="n", col=color.id[2]) #********* plot the fits of log(probability density) ************* # windows() par(pty="s", font=6, font.main=6, cex.lab=1.2, cex.main=1.4) plot(log.x.high.1, value.straight.line.den.theo.high.1, type="l", pch=21, lwd=2, cex=1, col=color.id[1], main="Linear fits of log (p) - Fig. 5", ylab="Log (p)", xlab=expression(Log (V[L])), xlim=c(-2,3), ylim=c(-5,1)) legend(.7,.7, name.series[1], lwd=2, bty="n", col=color.id[1]) grid(nx = NULL, ny = NULL, col = "lightgray", lty = "dotted", lwd = par("lwd")) lines(get(paste("log.x.high.",2,sep="")), get(paste("value.straight.line.den.theo.high.",2,sep="")), type="l", pch=21, lwd=2, cex=1, col=color.id[2]) legend(.7,.2, name.series[2], lwd=2, bty="n", col=color.id[2]) #********* plot the frequency density ************* # windows() par(pty="s", font=6, font.main=6, cex.lab=1.2, cex.main=1.4) matplot(linear.x.1, freq.linear.1, type="l", log="xy", pch=21, lwd=2, cex=1, lty=1, col=color.id[1], main="Comparison of frequency densities - Fig. 6", ylab=expression(paste("f, ( ",m^{-3},")",sep="")), xlab=expression(paste(V[L], ", ( ",m^{3},")",sep="")),xlim=c(10^-5,10^3), ylim=c(10^-4,10^5)) legend(.1,10^4.5, name.series[1], lwd=2, bty="n", lty=1, col=color.id[1]) grid(nx = NULL, ny = NULL, col = "lightgray", lty = "dotted",lwd = par("lwd")) lines(get(paste("linear.x.",2,sep="")), get(paste("freq.linear.",2,sep="")),type="l", pch=21,lwd=2,cex=1,col=color.id[2]) legend(.1,10^4, name.series[2], lwd=2, bty="n", col=color.id[2]) #********* plot the fit of frequency density ************* # windows() par(pty="s", font=6, font.main=6, cex.lab=1.2, cex.main=1.4) plot(log.x.high.1, value.straight.line.freq.theo.high.1, type="l", pch=21, lwd=2, cex=1, col=color.id[1], main="Linear fits of log (f) - Fig. 7", ylab="Log (f)", xlab=expression(paste("Log (",V[L],")", sep="")), xlim=c(-2,3), ylim=c(-4,4)) legend(.7,3.5, name.series[1], lwd=2, bty="n", col=color.id[1]) grid(nx = NULL, ny = NULL, col = "lightgray", lty = "dotted",lwd = par("lwd")) lines(get(paste("log.x.high.",2,sep="")), get(paste("value.straight.line.freq.theo.high.",2,sep="")),type="l", pch=21, lwd=2, cex=1, col=color.id[2]) legend(.7,3, name.series[2], lwd=2, bty="n", col=color.id[2]) for (volume.series in 1:2) { # windows() par(pty="s", font=6, font.main=6, cex.lab=1.2, cex.main=1.4) plot(get(paste("log.x.",volume.series,sep="")), get(paste("freq.log.",volume.series,sep="")), main=paste("Fit of log (f) - Fig. 8.",volume.series, sep=""), xlab=expression(paste("Log (",V[L],")", sep="")), ylab="Log (f)", type="l", col="blue", lwd=2) lines(get(paste("log.x.high.",volume.series,sep="")),get(paste("value.straight.line.freq.theo.high.", volume.series,sep="")), col="red", lwd=2, lty=2) grid(nx = NULL, ny = NULL, col = "lightgray", lty = "dotted", lwd = par("lwd"), equilogs = TRUE) beta.value<-round(get(paste("coef.fit.straight.line.freq.high.",volume.series,sep=""))[2],2) legend(-.3, 2.5, name.series[volume.series], col="red", cex=1, bty="n") text(0.2, 1.8, expression(beta )) text(.7, 1.8, paste(" = ", beta.value, sep="")) } dev.off() #********* write the probability density fits ************************* slope.series<-matrix(nrow = 2, ncol=2) for (coef.dataset in 1:2) { slope.series[coef.dataset,]<-(round(get(paste("coef.fit.straight.line.high.",coef.dataset,sep="")),2)) } write(slope.series[,2],file="slope.dat",append=TRUE,ncolumns=2,sep="\t") } #************** stop loop on events *************