src_Case_Study/Case_Study_graphical_representation.R

In alexdidkovskyi/RDDOOK3: Random Domain Decompositions for Object-Oriented Kriging

source('src_Case_Study/Case_Study_package.R')
### Graphical representation of results
#Graphical part
X <- as.vector(as.matrix(read.table(paste0(path_to_data,'X.txt'), header = T)))
Y <- as.vector(as.matrix(read.table(paste0(path_to_data,'Y.txt'), header = T)))
lenX <- length(X)
lenY <- length(Y)
M <- 100


ds.pred <- in.poly.matrix
fpred.d <- t(clr2density.mv(t(aggregated_results), t, t_step))
mean.field  <- t(mean.nr(t(fpred.d), t, t_step))
ds.pred[which(!is.na(in.poly.matrix))] <- mean.field

tiff(paste0(path_to_data, 'Images/Kriging-K-',K,'-new.tiff'), width=800, height=800)
par(cex=1.5)
image.plot(X, Y, ds.pred, col = rev(tim.colors(M)),
           main = paste('Ordinary Kriging - K=',K,sep=''),
           zlim = c(0.5,10), asp=1)
plot(shape, add=T)
dev.off()

#bootstrap variance
boot.map <- in.poly.matrix
boot.map[which(!is.na(in.poly.matrix))] <- bootstrap_variance
tiff(paste0(path_to_data,'Images/Bootstrap-Var-K-',K,'.tiff'), width=800, height=800)
par(cex = 1.5)
image.plot(X, Y, boot.map, col = (tim.colors(M)),
           main = paste('Bootstrap variance - K=',K,sep=''),
           zlim = c(0,13), asp = 1)
plot(shape, add = T)
dev.off()



#Kriging variance
krigv.map <- in.poly.matrix
krigv.map[which(!is.na(in.poly.matrix))] <-  kriging_variance
tiff(paste0(path_to_data,'Images/Kriging-Var-K-',K,'.tiff'),width=800, height=800)
par(cex = 1.5)
image.plot(X,Y,krigv.map,col = (tim.colors(M)),
           main = paste('Kriging variance - K=',K ,sep = ''),
           zlim = c(0, 40.07), asp = 1)
plot(shape, add = T)
dev.off()



mean.field.graph <- mean.field
fpred.d.graph <- fpred.d


col.scale <- rev(tim.colors(M))
breaks <- quantile(mean.field.graph, (0:M)/M, na.rm = TRUE)
color <- rep(0,length(mean.field.graph))
for(i in 1:length(mean.field.graph)){
  color[i] <- col.scale[min(which(mean.field.graph[i] <= (breaks)),na.rm = T)]
}
pdf(paste0(path_to_data, 'Images/Predictions-K-',K,'.pdf'), width = 16, height = 8)
par(mfrow = c(1, 2), cex = 1.25)
ds.pred <- in.poly.matrix
ds.pred[which(!is.na(in.poly.matrix))] <- mean.field.graph
image.plot(X, Y, ds.pred,
           col = rev(tim.colors(M)), breaks = breaks,
           main = paste('Predicted Mean DO - K=', K, sep = ''),
           zlim = c(1,8), asp = 1)
points(grid[seq(1, ngrid, length.out = 50), 1],grid[seq(1, ngrid, length.out = 50), 2], pch = 1,
       cex = .75, col = 1)
plot(shape, add = T)

par(cex = 1.25)
matplot(t, t(fpred.d.graph[seq(1, ngrid, length.out = 50), ]), ylim = c(0, 0.9),
        col = color[seq(1, ngrid, length.out = 50)], type='l', lty=1,
        ylab = 'Predicted densities', xlab = 'DO', main = paste('Predicted Densities of DO - K=',K,sep = ''))
dev.off()


#------- Plots of variances

  pdf(paste0(path_to_data,'Images/Variances-K-',K,'.pdf'), width = 16, height = 8)
  par(mfrow = c(1, 2), cex = 1.25, mar = c(5, 4, 4, 3) + 0.1)
  image.plot(X, Y, boot.map,col = (tim.colors(M)),
             main = paste('Bootstrap variance - K=', K, sep = ''),
             zlim = c(0, 13), asp = 1)
  plot(shape, add = T)

  par(cex = 1.25)
  image.plot(X, Y, krigv.map, col = (tim.colors(M)),
             main = paste('Kriging variance - K=', K, sep = ''),
             zlim = c(0, 41), asp = 1)
  plot(shape, add=T)
  dev.off()

#------- Plots of quartiles
# Select K (user)
K <- 8
N <- dim(fpred.d.graph-5)
quant.DO.1 = c(rep(NA,5), quantile.mv(t, fpred.d.graph[-(1:5),], .25))
quant.DO.2 = c(rep(NA,5), quantile.mv(t, fpred.d.graph[-(1:5),], .5))
quant.DO.3 = c(rep(NA,5), quantile.mv(t, fpred.d.graph[-(1:5),], .75))

quant1 <- in.poly.matrix
quant1[which(!is.na(in.poly.matrix))] <- quant.DO.1
quant2 <- in.poly.matrix
quant2[which(!is.na(in.poly.matrix))] <- quant.DO.2
quant3 <- in.poly.matrix
quant3[which(!is.na(in.poly.matrix))] <- quant.DO.3


pdf(paste0(path_to_data,'Images/Results-quantile',K,'.pdf'),width=16, height=8)
par(mfrow <- c(1,3), cex <- 1.5, mar <- c(5, 4, 4, 3) + 0.1)
image.plot(X, Y, quant1,col = rev(tim.colors(M)),
           main ='(a) 1st Predicted Quartile',
           zlim = c(0, 8.5), asp = 1, axes = F, xlab = '', ylab = '')
plot(shape, add = T)

par(cex = 1.5, mar = c(5, 4, 4, 3) + 0.1)
image.plot(X,Y,quant2,col = rev(tim.colors(M)),
           main = '(b) 2nd Predicted Quartile',
           zlim = c(0,8.5), asp = 1, axes = F, xlab = '', ylab = '')
plot(shape, add = T)

par(cex = 1.5, mar = c(5, 4, 4, 3) + 0.1)
image.plot(X, Y, quant3,col = rev(tim.colors(M)),
           main = '(c) 3rd Predicted Quartile',
           zlim = c(0, 8.5), asp = 1, axes = F, xlab = '', ylab = '')
plot(shape, add = T)
dev.off()


#------- Plots of probability of having DO < 2mg/l
# Set K (user)
K <- 8
N <- dim(fpred.d.graph - 5)
p.pred.2m <- c(rep(NA, 5), p.dens.mv(t, fpred.d.graph[-(1:5), ], 2))

pp1 <- in.poly.matrix
pp1[which(!is.na(in.poly.matrix))] <- p.pred.2m

pp.class <- pp1
pp.class[which(!is.na(in.poly.matrix))] <- ifelse(pp1[which(!is.na(in.poly.matrix))] > 0.5, 1, 0)

pdf(paste0(path_to_data,'Images/Ppred-cl-',K,'.pdf'), width = 14, height = 8)
par(mfrow = c(1,2), cex = 1.5, mar = c(5, 4, 4, 3.5) + 0.1)
image.plot(X, Y, pp1, col = (tim.colors(M)),
           main = paste0('(a) Probability of DO < 2 mg/l (K=', K,')'),
           zlim = c(0,1), asp = 1, axes = F, xlab = '', ylab = '')
plot(shape, add = T)

par(cex = 1.5, mar = c(5, 4, 4, 3.5) + 0.1)
image.plot(X, Y, pp.class,col = c('blue', 'red'),
           main = paste0('(b) Dead zones (K=',K,')'), breaks = c(0,0.5,1),
           zlim = c(0,1), asp = 1, axes = F, xlab = '', ylab = '', lab.breaks = c(0,0.5,1))
plot(shape, add = T)

dev.off()