From 7a6d75e7d6edbc6d540cc2bfaeb9b6a3a7bd8cf4 Mon Sep 17 00:00:00 2001 From: reedmaxwell Date: Fri, 7 Sep 2018 13:30:05 -0600 Subject: [PATCH] Update Waterfall_4panel_Plot.R --- .../R_scripts/Waterfall_4panel_Plot.R | 328 +++++++++--------- 1 file changed, 164 insertions(+), 164 deletions(-) diff --git a/Examples/Hillslope_Simulations/paper_cases/R_scripts/Waterfall_4panel_Plot.R b/Examples/Hillslope_Simulations/paper_cases/R_scripts/Waterfall_4panel_Plot.R index fb472de..909c811 100644 --- a/Examples/Hillslope_Simulations/paper_cases/R_scripts/Waterfall_4panel_Plot.R +++ b/Examples/Hillslope_Simulations/paper_cases/R_scripts/Waterfall_4panel_Plot.R @@ -1,167 +1,167 @@ -### Waterfall and ET Hillslope plots -### LEC, RMM 25-May-18, (edited 7-Sep-18) - -##This script loads in EcoSLIM output and makes Figures 4 and 5 from Maxwell et al Ecohydrology 2018 - -library("MASS") +### Waterfall and ET Hillslope plots +### LEC, RMM 25-May-18, (edited 7-Sep-18) + +##This script loads in EcoSLIM output and makes Figures 4 and 5 from Maxwell et al Ecohydrology 2018 + +library("MASS") library(fields) rm(list=ls()) - -## you may need to edit the working directory -setwd('~/EcoSLIM/Examples/Hillslope_Simulations/paper_cases/R_scripts') + +## you may need to edit the working directory uncomment and adjust as needed +##setwd('~/EcoSLIM/Examples/Hillslope_Simulations/paper_cases/R_scripts') source("./slicedens.R") - -months=c("OCT","NOV","DEC","JAN","FEB","MAR","APR","MAY","JUN","JUL","AUG","SEP") - -## make waterfall four panel (Fig 4) -fout = 'Figure4_waterfall_4-panel.pdf' -pdf(file=fout,width=9, height=9) - -par(mfrow=c(2,2)) - -#Read EcoSlim output files and put in independent matrices for plotting -## load trees with ER forcing -fin='./ER_hillslope_trees/SLIM_hillslope_ER_trees_exited_particles.bin' -length=file.info(fin)$size/8 -ncol=8 -print(length)/ncol -to.read = file(fin,'rb') -part=matrix(0, nrow=(length/ncol), ncol=ncol) -for(i in 1:(length/ncol)){ - part[i,]= readBin(to.read, double(), endian='little', size=8,n=ncol) -} -close(to.read) -colnames(part)=c("Time", "x", "y", "z", "Ptime", "Mass", "Comp", "ExitStatus") -npart=nrow(part) - -#3D Histograms of particle ages -out = which( (part[,8] == 1) & (part[,1]>=(19*8760))) -et = which((part[,8] == 2) & (part[,1]>=(19*8760))) -den3d_out_er_trees <- kde2d( (part[out,5])/24/365,part[out,1]/24/365, n = 50) - -den3d_et_er_trees_time <- kde2d( (part[et,5])/24/365,part[et,1]/24/365, n = 50) -den3d_et_er_trees_space <- kde2d( (part[et,5])/24/365,part[et,2], n = 100) - -#Outflow -x=part[out,5]/24/365 -y=part[out,1]/24/365 -z <- y -trans=0.7 -fcol <- rbind(c(0,.1,.5,trans), c(.3,.8,.8,trans), c(1,1,0, trans)) -fcol=fcol[3:1,] -lcol <- rbind(c(0,.3,.3,.8), c(.1,.1,.2,.7), c(0,0,1,.65)) - - -slicedens(x,y,z,fcol=fcol, bcol='white', lcol=lcol,gboost=.9, yinc=1, slices=12, xlab="Outflow Residence Time [YR] ",ylab="",main="A) ER Trees Outflow",yaxlab=months) - -#ET -x=part[et,5]/24/365 -y=part[et,1]/24/365 -z <- y - -slicedens(x,y,z,fcol=fcol, bcol='white', lcol=lcol,gboost=.90, yinc=7, slices=12, xlab="ET Residence Time [YR]",ylab="", main="B) ER Trees ET",yaxlab=months) - - -## load shrubs with LW forcing -fin='../LW_hillslope_shrub/SLIM_hillslope_LW_shrub_exited_particles.bin' -length=file.info(fin)$size/8 -ncol=8 -print(length)/ncol -to.read = file(fin,'rb') -part=matrix(0, nrow=(length/ncol), ncol=ncol) -for(i in 1:(length/ncol)){ - part[i,]= readBin(to.read, double(), endian='little', size=8,n=ncol) -} -close(to.read) -colnames(part)=c("Time", "x", "y", "z", "Ptime", "Mass", "Comp", "ExitStatus") -npart=nrow(part) - -#3D Histograms of particle ages -out = which( (part[,8] == 1) & (part[,1]>=(19*8760))) -et = which((part[,8] == 2) & (part[,1]>=(19*8760))) -den3d_out_lw_shrubs <- kde2d( (part[out,5])/24/365,part[out,1]/24/365, n = 50) -den3d_et_lw_shrubs_time <- kde2d( (part[et,5])/24/365,part[et,1]/24/365, n = 50) -den3d_et_lw_shrubs_space <- kde2d( (part[et,5])/24/365,part[et,2], n = 100) - - -## load trees with LW forcing -fin='../LW_hillslope_trees/SLIM_hillslope_LW_trees_exited_particles.bin' -length=file.info(fin)$size/8 -ncol=8 -print(length)/ncol -to.read = file(fin,'rb') -part=matrix(0, nrow=(length/ncol), ncol=ncol) -for(i in 1:(length/ncol)){ - part[i,]= readBin(to.read, double(), endian='little', size=8,n=ncol) -} -close(to.read) -colnames(part)=c("Time", "x", "y", "z", "Ptime", "Mass", "Comp", "ExitStatus") -npart=nrow(part) - -#3D Histograms of particle ages -out = which( (part[,8] == 1) & (part[,1]>=(19*8760))) -et = which((part[,8] == 2) & (part[,1]>=(19*8760))) -den3d_out_lw_trees <- kde2d( (part[out,5])/24/365,part[out,1]/24/365, n = 50) -den3d_et_lw_trees_time <- kde2d( (part[et,5])/24/365,part[et,1]/24/365, n = 50) -den3d_et_lw_trees_space <- kde2d( (part[et,5])/24/365,part[et,2], n = 100) - -#Outflow -x=part[out,5]/24/365 -y=part[out,1]/24/365 -z <- y -trans=0.7 -fcol <- rbind(c(0,.1,.5,trans), c(.3,.8,.8,trans), c(1,1,0, trans)) -fcol=fcol[3:1,] -lcol <- rbind(c(0,.3,.3,.8), c(.1,.1,.2,.7), c(0,0,1,.65)) - - -slicedens(x,y,z,fcol=fcol, bcol='white', lcol=lcol,gboost=.9, yinc=1, slices=12, xlab="Outflow Residence Time [YR]",main="C) LW Trees Outflow",ylab="",yaxlab=months) - -#ET -x=part[et,5]/24/365 -y=part[et,1]/24/365 -z <- y - -slicedens(x,y,z,fcol=fcol, bcol='white', lcol=lcol,gboost=.90, yinc=7, slices=12, xlab="ET Residence Time [YR]",main="D) LW Trees ET",ylab="",yaxlab=months) - -dev.off() - -## load shrubs with ER forcing -fin='./ER_hillslope_shrub/SLIM_hillslope_ER_shrub_exited_particles.bin' -length=file.info(fin)$size/8 -ncol=8 -print(length)/ncol -to.read = file(fin,'rb') -part=matrix(0, nrow=(length/ncol), ncol=ncol) -for(i in 1:(length/ncol)){ - part[i,]= readBin(to.read, double(), endian='little', size=8,n=ncol) -} -close(to.read) -colnames(part)=c("Time", "x", "y", "z", "Ptime", "Mass", "Comp", "ExitStatus") -npart=nrow(part) - -#3D Histograms of particle ages -out = which( (part[,8] == 1) & (part[,1]>=(19*8760))) -et = which((part[,8] == 2) & (part[,1]>=(19*8760))) -den3d_out_er_shrubs <- kde2d( (part[out,5])/24/365,part[out,1]/24/365, n = 50) -den3d_et_er_shrubs_time <- kde2d( (part[et,5])/24/365,part[et,1]/24/365, n = 50) -den3d_et_er_shrubs_space <- kde2d( (part[et,5])/24/365,part[et,2], n = 100) - -## make ET distribution four panel (Fig 5) -fout = 'Figure5_ET_hillsllope_4-panel.pdf' -pdf(file=fout) - -par(mfrow=c(2,2)) - -image(den3d_et_er_shrubs_space,xlab = "ET residence time [yr]", ylab = "X Location Hillslope [m]",col = heat.colors(9),zlim=c(0.0001,0.05),xlim=c(0,10),breaks=c(0.000001,0.000005,0.00001,0.00005,0.0001,0.0005,0.001,0.005,0.01,0.05),main="A) ER Shrubs") - -image(den3d_et_er_trees_space,xlab = "ET residence time [yr]", ylab = "X Location Hillslope [m]",col = heat.colors(9),zlim=c(0.0001,0.05),xlim=c(0,10),breaks=c(0.000001,0.000005,0.00001,0.00005,0.0001,0.0005,0.001,0.005,0.01,0.05),main="B) ER Trees") - -image(den3d_et_lw_shrubs_space,xlab = "ET residence time [yr]", ylab = "X Location Hillslope [m]",col = heat.colors(9),zlim=c(0.0001,0.05),xlim=c(0,10),breaks=c(0.000001,0.000005,0.00001,0.00005,0.0001,0.0005,0.001,0.005,0.01,0.05),main="C) LW Shrubs") - -image(den3d_et_lw_trees_space,xlab = "ET residence time [yr]", ylab = "X Location Hillslope [m]",col = heat.colors(9),zlim=c(0.0001,0.05),xlim=c(0,10),breaks=c(0.000001,0.000005,0.00001,0.00005,0.0001,0.0005,0.001,0.005,0.01,0.05),main="D) LW Trees") - -image.plot(den3d_et_lw_trees_space, zlim=c(0.0001,0.05), cex=0.7,col=(heat.colors(9)), legend.only=T, legend.width=0.5, legend.shrink=0.5, legend.args=list(text=' Density [-]',cex=0.85,side=3,line=0.1), smallplot=c(.65,.68, .55,.75)) -dev.off() - - + +months=c("OCT","NOV","DEC","JAN","FEB","MAR","APR","MAY","JUN","JUL","AUG","SEP") + +## make waterfall four panel (Fig 4) +fout = 'Figure4_waterfall_4-panel.pdf' +pdf(file=fout,width=9, height=9) + +par(mfrow=c(2,2)) + +#Read EcoSlim output files and put in independent matrices for plotting +## load trees with ER forcing +fin='./ER_hillslope_trees/SLIM_hillslope_ER_trees_exited_particles.bin' +length=file.info(fin)$size/8 +ncol=8 +print(length)/ncol +to.read = file(fin,'rb') +part=matrix(0, nrow=(length/ncol), ncol=ncol) +for(i in 1:(length/ncol)){ + part[i,]= readBin(to.read, double(), endian='little', size=8,n=ncol) +} +close(to.read) +colnames(part)=c("Time", "x", "y", "z", "Ptime", "Mass", "Comp", "ExitStatus") +npart=nrow(part) + +#3D Histograms of particle ages +out = which( (part[,8] == 1) & (part[,1]>=(19*8760))) +et = which((part[,8] == 2) & (part[,1]>=(19*8760))) +den3d_out_er_trees <- kde2d( (part[out,5])/24/365,part[out,1]/24/365, n = 50) + +den3d_et_er_trees_time <- kde2d( (part[et,5])/24/365,part[et,1]/24/365, n = 50) +den3d_et_er_trees_space <- kde2d( (part[et,5])/24/365,part[et,2], n = 100) + +#Outflow +x=part[out,5]/24/365 +y=part[out,1]/24/365 +z <- y +trans=0.7 +fcol <- rbind(c(0,.1,.5,trans), c(.3,.8,.8,trans), c(1,1,0, trans)) +fcol=fcol[3:1,] +lcol <- rbind(c(0,.3,.3,.8), c(.1,.1,.2,.7), c(0,0,1,.65)) + + +slicedens(x,y,z,fcol=fcol, bcol='white', lcol=lcol,gboost=.9, yinc=1, slices=12, xlab="Outflow Residence Time [YR] ",ylab="",main="A) ER Trees Outflow",yaxlab=months) + +#ET +x=part[et,5]/24/365 +y=part[et,1]/24/365 +z <- y + +slicedens(x,y,z,fcol=fcol, bcol='white', lcol=lcol,gboost=.90, yinc=7, slices=12, xlab="ET Residence Time [YR]",ylab="", main="B) ER Trees ET",yaxlab=months) + + +## load shrubs with LW forcing +fin='../LW_hillslope_shrub/SLIM_hillslope_LW_shrub_exited_particles.bin' +length=file.info(fin)$size/8 +ncol=8 +print(length)/ncol +to.read = file(fin,'rb') +part=matrix(0, nrow=(length/ncol), ncol=ncol) +for(i in 1:(length/ncol)){ + part[i,]= readBin(to.read, double(), endian='little', size=8,n=ncol) +} +close(to.read) +colnames(part)=c("Time", "x", "y", "z", "Ptime", "Mass", "Comp", "ExitStatus") +npart=nrow(part) + +#3D Histograms of particle ages +out = which( (part[,8] == 1) & (part[,1]>=(19*8760))) +et = which((part[,8] == 2) & (part[,1]>=(19*8760))) +den3d_out_lw_shrubs <- kde2d( (part[out,5])/24/365,part[out,1]/24/365, n = 50) +den3d_et_lw_shrubs_time <- kde2d( (part[et,5])/24/365,part[et,1]/24/365, n = 50) +den3d_et_lw_shrubs_space <- kde2d( (part[et,5])/24/365,part[et,2], n = 100) + + +## load trees with LW forcing +fin='../LW_hillslope_trees/SLIM_hillslope_LW_trees_exited_particles.bin' +length=file.info(fin)$size/8 +ncol=8 +print(length)/ncol +to.read = file(fin,'rb') +part=matrix(0, nrow=(length/ncol), ncol=ncol) +for(i in 1:(length/ncol)){ + part[i,]= readBin(to.read, double(), endian='little', size=8,n=ncol) +} +close(to.read) +colnames(part)=c("Time", "x", "y", "z", "Ptime", "Mass", "Comp", "ExitStatus") +npart=nrow(part) + +#3D Histograms of particle ages +out = which( (part[,8] == 1) & (part[,1]>=(19*8760))) +et = which((part[,8] == 2) & (part[,1]>=(19*8760))) +den3d_out_lw_trees <- kde2d( (part[out,5])/24/365,part[out,1]/24/365, n = 50) +den3d_et_lw_trees_time <- kde2d( (part[et,5])/24/365,part[et,1]/24/365, n = 50) +den3d_et_lw_trees_space <- kde2d( (part[et,5])/24/365,part[et,2], n = 100) + +#Outflow +x=part[out,5]/24/365 +y=part[out,1]/24/365 +z <- y +trans=0.7 +fcol <- rbind(c(0,.1,.5,trans), c(.3,.8,.8,trans), c(1,1,0, trans)) +fcol=fcol[3:1,] +lcol <- rbind(c(0,.3,.3,.8), c(.1,.1,.2,.7), c(0,0,1,.65)) + + +slicedens(x,y,z,fcol=fcol, bcol='white', lcol=lcol,gboost=.9, yinc=1, slices=12, xlab="Outflow Residence Time [YR]",main="C) LW Trees Outflow",ylab="",yaxlab=months) + +#ET +x=part[et,5]/24/365 +y=part[et,1]/24/365 +z <- y + +slicedens(x,y,z,fcol=fcol, bcol='white', lcol=lcol,gboost=.90, yinc=7, slices=12, xlab="ET Residence Time [YR]",main="D) LW Trees ET",ylab="",yaxlab=months) + +dev.off() + +## load shrubs with ER forcing +fin='./ER_hillslope_shrub/SLIM_hillslope_ER_shrub_exited_particles.bin' +length=file.info(fin)$size/8 +ncol=8 +print(length)/ncol +to.read = file(fin,'rb') +part=matrix(0, nrow=(length/ncol), ncol=ncol) +for(i in 1:(length/ncol)){ + part[i,]= readBin(to.read, double(), endian='little', size=8,n=ncol) +} +close(to.read) +colnames(part)=c("Time", "x", "y", "z", "Ptime", "Mass", "Comp", "ExitStatus") +npart=nrow(part) + +#3D Histograms of particle ages +out = which( (part[,8] == 1) & (part[,1]>=(19*8760))) +et = which((part[,8] == 2) & (part[,1]>=(19*8760))) +den3d_out_er_shrubs <- kde2d( (part[out,5])/24/365,part[out,1]/24/365, n = 50) +den3d_et_er_shrubs_time <- kde2d( (part[et,5])/24/365,part[et,1]/24/365, n = 50) +den3d_et_er_shrubs_space <- kde2d( (part[et,5])/24/365,part[et,2], n = 100) + +## make ET distribution four panel (Fig 5) +fout = 'Figure5_ET_hillsllope_4-panel.pdf' +pdf(file=fout) + +par(mfrow=c(2,2)) + +image(den3d_et_er_shrubs_space,xlab = "ET residence time [yr]", ylab = "X Location Hillslope [m]",col = heat.colors(9),zlim=c(0.0001,0.05),xlim=c(0,10),breaks=c(0.000001,0.000005,0.00001,0.00005,0.0001,0.0005,0.001,0.005,0.01,0.05),main="A) ER Shrubs") + +image(den3d_et_er_trees_space,xlab = "ET residence time [yr]", ylab = "X Location Hillslope [m]",col = heat.colors(9),zlim=c(0.0001,0.05),xlim=c(0,10),breaks=c(0.000001,0.000005,0.00001,0.00005,0.0001,0.0005,0.001,0.005,0.01,0.05),main="B) ER Trees") + +image(den3d_et_lw_shrubs_space,xlab = "ET residence time [yr]", ylab = "X Location Hillslope [m]",col = heat.colors(9),zlim=c(0.0001,0.05),xlim=c(0,10),breaks=c(0.000001,0.000005,0.00001,0.00005,0.0001,0.0005,0.001,0.005,0.01,0.05),main="C) LW Shrubs") + +image(den3d_et_lw_trees_space,xlab = "ET residence time [yr]", ylab = "X Location Hillslope [m]",col = heat.colors(9),zlim=c(0.0001,0.05),xlim=c(0,10),breaks=c(0.000001,0.000005,0.00001,0.00005,0.0001,0.0005,0.001,0.005,0.01,0.05),main="D) LW Trees") + +image.plot(den3d_et_lw_trees_space, zlim=c(0.0001,0.05), cex=0.7,col=(heat.colors(9)), legend.only=T, legend.width=0.5, legend.shrink=0.5, legend.args=list(text=' Density [-]',cex=0.85,side=3,line=0.1), smallplot=c(.65,.68, .55,.75)) +dev.off() + +