Rutils/temp/plot_yearly2.r
In manfredo89/ED2io: Manages Input/Output for Ecosystem Demography 2 Model
#==========================================================================================#
#==========================================================================================#
# Leave these commands at the beginning. They will refresh the session. #
#------------------------------------------------------------------------------------------#
rm(list=ls())
graphics.off()
#==========================================================================================#
#==========================================================================================#
#==========================================================================================#
#==========================================================================================#
# Here is the user defined variable section. #
#------------------------------------------------------------------------------------------#
#----- Paths. -----------------------------------------------------------------------------#
args <- commandArgs(TRUE)
arg.runtype <- as.character(args[])
if (length(args) > 1){
cat("Only one directory is supported now \n")
stop()
}
if (length(args)==0) {
cat("No arguments were passed, defaulting to lianas \n")
arg.runtype="lianas"
}
run.type = arg.runtype
here = getwd() # Current directory.
there = paste(here,"/../",run.type[1],sep='') # Directory where analyses/history are
srcdir = here # Source directory.
outroot = paste(there,"/figures",sep='') # Directory for figures
cat("here:",here,"\n")
cat("there:",there,"\n")
cat("srcdir:",srcdir,"\n")
cat("outroot:",outroot,"\n")
#------------------------------------------------------------------------------------------#
#----- Time options. ----------------------------------------------------------------------#
monthbeg = 01 # First month to use
yearbeg = 2000 # First year to consider
yearend = 2050 # Maximum year to consider
reload.data = T # Should I reload partially loaded data?
#sasmonth.short = c(2,5,8,11) # Months for SAS plots (short runs)
sasmonth.short = c(2,4,6,8,10,12) # Months for SAS plots (short runs)
sasmonth.long = c(2,4,6,8,10,12) # Months for SAS plots (short runs)
#sasmonth.long = 5 # Months for SAS plots (long runs)
#nyears.long = 800 # Runs longer than this are considered long runs.
#------------------------------------------------------------------------------------------#
#----- Name of the simulations. -----------------------------------------------------------#
myplaces = c("paracou")
#------------------------------------------------------------------------------------------#
#----- Plot options. ----------------------------------------------------------------------#
outform = c("pdf") # Formats for output file. Supported formats are:
# - "X11" - for printing on screen
# - "quartz" - for printing on Mac OS screen
# - "eps" - for postscript printing
# - "png" - for PNG printing
# - "tif" - for TIFF printing
# - "pdf" - for PDF printing
depth = 96 # PNG resolution, in pixels per inch
paper = "letter" # Paper size, to define the plot shape
ptsz = 16 # Font size.
lwidth = 2.5 # Line width
plotgrid = TRUE # Should I plot the grid in the background?
sasfixlimits = FALSE # Use a fixed scale for size and age-structure
# plots? (FALSE will set a suitable scale for
# each plot)
fcgrid = TRUE # Include a grid on the filled contour plots?
ncolshov = 200 # Target number of colours for Hovmoller diagrams.
hovgrid = TRUE # Include a grid on the Hovmoller plots?
legwhere = "topleft" # Where should I place the legend?
inset = 0.01 # Inset between legend and edge of plot region.
scalleg = 0.40 # Expand y limits by this relative amount to fit
# the legend
cex.main = 0.8 # Scale coefficient for the title
theta = 315. # Azimuth for perspective projection
phi = 30. # Vertical angle for perspective projection
ltheta = -210. # Azimuth angle for light
shade = 0.125 # Shade intensity
expz = 0.5 # Expansion factor for Z axis
cexmin = 0.5 # Minimum "head" size of the lollipop
cexmax = 3.0 # Maximum "head" size of the lollipop
ylnudge = 0.05 # Nudging factor for ylimit
ptype = "l" # Type of plot
ptyped = "p" # Type of plot
ptypeb = "o" # Type of plot
drought.mark = FALSE # Put a background to highlight droughts?
drought.yeara = 1605 # First year that has drought
drought.yearz = 1609 # Last year that has drought
months.drought = c(12,1,2,3) # Months with drought
ibackground = 0 # Background settings (check load_everything.r)
plot.nplant.hystogram = T # Plot custom dbh distribution (works but bad implementation, correct or remove)
options(bitmapType='cairo')
#------------------------------------------------------------------------------------------#
#------ Miscellaneous settings. -----------------------------------------------------------#
slz.min = -5.0 # The deepest depth that trees access water.
idbh.type = 3 # Type of DBH class
# 1 -- Every 10 cm until 100cm; > 100cm
# 2 -- 0-10; 10-20; 20-35; 35-50; 50-70; > 70 (cm)
# 3 -- 0-10; 10-35; 35-55; > 55 (cm)
klight = 0.8 # Weighting factor for maximum carbon balance
corr.growth.storage = 1.0 # Correction factor to be applied to growth and
# storage respiration
iallom = 3 # Allometry to use
#------------------------------------------------------------------------------------------#
#==========================================================================================#
#==========================================================================================#
#==========================================================================================#
#==========================================================================================#
#==========================================================================================#
#==========================================================================================#
# NO NEED TO CHANGE ANYTHING BEYOND THIS POINT UNLESS YOU ARE DEVELOPING THE CODE... #
#==========================================================================================#
#==========================================================================================#
#==========================================================================================#
#==========================================================================================#
#==========================================================================================#
#==========================================================================================#
analy.path = paste(there,"analy",sep='/') # analysis folder
com.list = list.files(analy.path, pattern = "*-Q-*")
if(length(run.type) > 1) com.list = com.list[duplicated(com.list)]
if(length(com.list) < 24){ stop("ERROR: two few files in analysis folder") }
#------------------------------------------------------------------------------------------#
# Keep only full years #
#------------------------------------------------------------------------------------------#
# first year of available data in the folder
if (is.na(yearbeg)) yearbeg = min(sub('\\-.*','',sub("^.*?Q-","",com.list)))
# last year if available data in the folder
if (is.na(yearend)) yearend = max(sub('\\-.*','',sub("^.*?Q-","",com.list)))
yearbeg = as.integer(yearbeg) # convert yeara to factor
yearend = as.integer(yearend) # convert yearz to factor
#----- Load some packages and scripts. ----------------------------------------------------#
source(file.path(srcdir,"load.everything.r"))
#------------------------------------------------------------------------------------------#
#----- Set how many formats we must output. -----------------------------------------------#
outform = tolower(outform)
nout = length (outform)
#------------------------------------------------------------------------------------------#
#----- Avoid unecessary and extremely annoying beeps. -------------------------------------#
options(locatorBell=FALSE)
#------------------------------------------------------------------------------------------#
#----- Load observations. -----------------------------------------------------------------#
#obsrfile = file.path(srcdir,"LBA_MIP.v8.RData")
#load(file=obsrfile)
#----- Define plot window size ------------------------------------------------------------#
size = plotsize(proje=FALSE,paper=paper)
#------------------------------------------------------------------------------------------#
#---- Create the main output directory in case there is none. -----------------------------#
if (! file.exists(outroot)) dir.create(outroot)
#------------------------------------------------------------------------------------------#
#------------------------------------------------------------------------------------------#
# Big place loop starts here... #
#------------------------------------------------------------------------------------------#
#for (place in myplaces){
place="paracou"
#----- Retrieve default information about this place and set up some variables. --------#
thispoi = locations(where=place,here=there,yearbeg=yearbeg,yearend=yearend
,monthbeg=monthbeg)
inpref = thispoi$pathin
outmain = file.path(outroot,place)
outpref = file.path(outmain,"yearly")
lieu = thispoi$lieu
iata = thispoi$iata
suffix = thispoi$iata
yeara = thispoi$yeara
yearz = thispoi$yearz
meszz = thispoi$monz
#---------------------------------------------------------------------------------------#
# Make sure we only deal with full years. #
#---------------------------------------------------------------------------------------#
if (monthbeg > 1) yeara = yeara + 1
if (meszz < 12) yearz = yearz - 1
monthbeg = 1
meszz = 12
if (yeara > yearz){
cat(" - Yeara: ",yeara,"\n")
cat(" - Yearz: ",yearz,"\n")
cat(" - Prefix: ",inpref,"\n")
cat(" - Invalid years, will not process data...","\n")
q("no")
}#end if
#---------------------------------------------------------------------------------------#
#----- Create the directories in case they don't exist. --------------------------------#
if (! file.exists(outmain)) dir.create(outmain)
if (! file.exists(outpref)) dir.create(outpref)
#---------------------------------------------------------------------------------------#
#----- Decide how frequently the cohort-level variables should be saved. ---------------#
if ((yearend - yearbeg + 1) <= nyears.long){
sasmonth = sasmonth.short
plot.ycomp = TRUE
}else{
sasmonth = sasmonth.long
# plot.ycomp = FALSE
plot.ycomp = TRUE
}#end if
#---------------------------------------------------------------------------------------#
#----- Print a banner to entretain the user. -------------------------------------------#
cat(" + Post-processing output from ",lieu,"...","\n")
#---------------------------------------------------------------------------------------#
# Flush all variables that will hold the data. #
#---------------------------------------------------------------------------------------#
ntimes = (yearz-yeara-1)*12+meszz+(12-monthbeg+1)
nyears = yearz-yeara+1
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#
# Make the RData file name, then we check whether we must read the files again #
# or use the stored RData. Notice that the path is the same for plot_ycomp.r and #
# plot_monthly, so you don't need to read in the data twice. #
#---------------------------------------------------------------------------------------#
path.data = file.path(there,"rdata_month")
if (! file.exists(path.data)) dir.create(path.data)
ed22.rdata = file.path(path.data,paste("M",place,"RData",sep="."))
ed22.status = file.path(path.data,paste("Mstatus_",place,".txt",sep=""))
if (reload.data && file.exists(ed22.rdata)){
#----- Load the modelled dataset. ---------------------------------------------------#
cat(" - Loading previous session...","\n")
load(ed22.rdata)
tresume = datum$ntimes + 1
if (ntimes > datum$ntimes){
datum = update.monthly( new.ntimes = ntimes
, old.datum = datum
, montha = monthbeg
, yeara = yeara
, inpref = inpref
, slz.min = slz.min
)#end update.monthly
}#end if
#------------------------------------------------------------------------------------#
}else{
cat(" - Starting new session...","\n")
tresume = 1
datum = create.monthly( ntimes = ntimes
, montha = monthbeg
, yeara = yeara
, inpref = inpref
, slz.min = slz.min
)#end create.monthly
}#end if
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#
# Check whether we have anything to update. #
#---------------------------------------------------------------------------------------#
complete = tresume > ntimes
#---------------------------------------------------------------------------------------#
#----- Copy some dimensions to scalars. ------------------------------------------------#
nzg = datum$nzg
nzs = datum$nzs
ndcycle = datum$ndcycle
isoilflg = datum$isoilflg
slz = datum$slz
slxsand = datum$slxsand
slxclay = datum$slxclay
ntext = datum$ntext
soil.prop = datum$soil.prop
dslz = datum$dslz
soil.depth = datum$soil.depth
soil.dry = datum$soil.dry
soil.poro = datum$soil.poro
ka = datum$ka
kz = datum$kz
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#
# Loop over all times in case there is anything new to be read. #
#---------------------------------------------------------------------------------------#
if (! complete){
#------------------------------------------------------------------------------------#
# This function will read the files. #
#------------------------------------------------------------------------------------#
datum = read.q.files(datum=datum,ntimes=ntimes,tresume=tresume,sasmonth=sasmonth)
#------------------------------------------------------------------------------------#
#------ Save the data to the R object. ----------------------------------------------#
cat(" + Saving data to ",basename(ed22.rdata),"...","\n")
save(datum,file=ed22.rdata)
#------------------------------------------------------------------------------------#
}#end if (! complete)
#---------------------------------------------------------------------------------------#
#----- Update status file with latest data converted into R. ---------------------------#
#latest = paste(datum$year[ntimes],datum$month[ntimes],sep=" ")
#dummy = write(x=latest,file=ed22.status,append=FALSE)
#---------------------------------------------------------------------------------------#
#----- Make some shorter versions of some variables. -----------------------------------#
mfac = datum$month
yfac = datum$year
emean = datum$emean
emsqu = datum$emsqu
qmean = datum$qmean
qmsqu = datum$qmsqu
szpft = datum$szpft
lu = datum$lu
patch = datum$patch
cohort = datum$cohort
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#
# Consolidate the yearly means for the long-term dynamics (the PFT and DBH/PFT #
# stuff). #
#---------------------------------------------------------------------------------------#
cat (" - Finding the annual statistics for multi-dimensional variables...","\n")
cat (" * Aggregating the annual mean of PFT-DBH variables...","\n")
for (vname in names(szpft)){
szpft[[vname]] = qapply(X=szpft[[vname]],INDEX=yfac,DIM=1,FUN=mean,na.rm=TRUE)
}#end for
#----- LU arrays. The "+1" column contains the total. --------------------------------#
cat (" * Aggregating the annual mean of LU variables...","\n")
for (vname in names(lu)){
lu [[vname]] = qapply(X=lu [[vname]],INDEX=yfac,DIM=1,FUN=mean,na.rm=TRUE)
}#end for
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#
# Here we find the monthly means for month, then compute the standard deviation. #
#---------------------------------------------------------------------------------------#
cat (" - Finding the monthly and annual means...","\n")
cat (" * Aggregating the monthly mean and standard deviation...","\n")
mmean = list()
msdev = list()
ymean = list()
ysdev = list()
for (vname in names(emean)){
if (vname %in% c("soil.temp","soil.water","soil.mstpot","soil.extracted")){
mmean[[vname]] = qapply(X=emean[[vname]], INDEX=mfac, DIM=1, FUN=mean, na.rm=TRUE)
msdev[[vname]] = qapply(X=emean[[vname]], INDEX=mfac, DIM=1, FUN=sd , na.rm=TRUE)
ymean[[vname]] = qapply(X=emean[[vname]], INDEX=yfac, DIM=1, FUN=mean, na.rm=TRUE)
ysdev[[vname]] = qapply(X=emean[[vname]], INDEX=yfac, DIM=1, FUN=sd , na.rm=TRUE)
}else if (vname %in% c("rain","runoff","intercepted","wshed")){
mmean[[vname]] = tapply(X=emean[[vname]], INDEX=mfac, FUN=mean, na.rm=TRUE)
msdev[[vname]] = tapply(X=emean[[vname]], INDEX=mfac, FUN=sd , na.rm=TRUE)
ymean[[vname]] = tapply(X=emean[[vname]], INDEX=yfac, FUN=sum , na.rm=TRUE)
ysdev[[vname]] = tapply(X=emean[[vname]], INDEX=yfac, FUN=sd , na.rm=TRUE)
}else{
mmean[[vname]] = tapply(X=emean[[vname]], INDEX=mfac, FUN=mean, na.rm=TRUE)
msdev[[vname]] = tapply(X=emean[[vname]], INDEX=mfac, FUN=sd , na.rm=TRUE)
ymean[[vname]] = tapply(X=emean[[vname]], INDEX=yfac, FUN=mean, na.rm=TRUE)
ysdev[[vname]] = tapply(X=emean[[vname]], INDEX=yfac, FUN=sd , na.rm=TRUE)
}#end if
#------------------------------------------------------------------------------------#
#----- Fix the bad data. ------------------------------------------------------------#
bad.mmean = ! is.finite(mmean[[vname]])
bad.msdev = ! is.finite(msdev[[vname]])
bad.ymean = ! is.finite(ymean[[vname]])
bad.ysdev = ! is.finite(ysdev[[vname]])
mmean[[vname]][bad.mmean] = NA
msdev[[vname]][bad.msdev] = 0.
ymean[[vname]][bad.ymean] = NA
ysdev[[vname]][bad.ysdev] = 0.
#------------------------------------------------------------------------------------#
}#end for
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#
# Here we find the Mean diurnal cycle for each month, then compute the standard #
# deviation. #
#---------------------------------------------------------------------------------------#
cat (" - Aggregating the annual mean and std. dev. of the diurnal cycle...","\n")
umean = list()
usdev = list()
for (vname in names(qmean)){
umean[[vname]] = qapply(qmean[[vname]],INDEX=yfac,DIM=1,FUN=mean,na.rm=TRUE)
usdev[[vname]] = qapply(qmean[[vname]],INDEX=yfac,DIM=1,FUN=sd ,na.rm=TRUE)
bad.umean = ! is.finite(umean[[vname]])
bad.usdev = ! is.finite(usdev[[vname]])
umean[[vname]][bad.umean] = NA
usdev[[vname]][bad.usdev] = 0.
}#end for
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#
# Remove all elements of the DBH/PFT class that do not have a single valid cohort #
# at any given time. #
#---------------------------------------------------------------------------------------#
empty = is.na(szpft$nplant) | szpft$nplant == 0
for (vname in names(szpft)) szpft[[vname]][empty] = NA
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#
# Convert mortality and recruitment so it is scaled between 0 and 100%. #
#---------------------------------------------------------------------------------------#
szpft$mort = 100. * (1.0 - exp(- szpft$mort ) )
szpft$dimort = 100. * (1.0 - exp(- szpft$dimort ) )
szpft$ncbmort = 100. * (1.0 - exp(- szpft$ncbmort ) )
szpft$recrpft = 100. * ( exp( szpft$recr ) - 1.0)
szpft$agb.mort = 100. * (1.0 - exp(- szpft$agb.mort ) )
szpft$agb.dimort = 100. * (1.0 - exp(- szpft$agb.dimort ) )
szpft$agb.ncbmort = 100. * (1.0 - exp(- szpft$agb.ncbmort ) )
szpft$agb.recrpft = 100. * ( exp( szpft$agb.recr ) - 1.0)
szpft$bsa.mort = 100. * (1.0 - exp(- szpft$bsa.mort ) )
szpft$bsa.dimort = 100. * (1.0 - exp(- szpft$bsa.dimort ) )
szpft$bsa.ncbmort = 100. * (1.0 - exp(- szpft$bsa.ncbmort ) )
szpft$bsa.recrpft = 100. * ( exp( szpft$bsa.recr ) - 1.0)
#---------------------------------------------------------------------------------------#
#----- Find which PFTs, land uses and transitions we need to consider ------------------#
pftave = apply( X = szpft$agb[,ndbh+1,]
, MARGIN = 2
, FUN = mean
, na.rm = TRUE
)#end apply
luave = apply( X = lu$agb
, MARGIN = 2
, FUN = mean
, na.rm = TRUE
)#end apply
distave = apply(X=lu$dist,MARGIN=c(2,3),FUN=mean)
selpft = is.finite(pftave ) & pftave > 0.
sellu = is.finite(luave ) & luave > 0.
seldist = is.finite(distave) & distave > 0.
n.selpft = sum(selpft )
n.sellu = sum(sellu )
n.seldist = sum(seldist)
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#
# Define a suitable scale for diurnal cycle... #
#---------------------------------------------------------------------------------------#
thisday = seq(from=0,to=ndcycle,by=1) * 24 / ndcycle
uplot = list()
uplot$levels = c(0,4,8,12,16,20,24)
uplot$n = 7
uplot$scale = "hours"
uplot$padj = rep(0,times=uplot$n)
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#
# Define a suitable scale for soil profile layers... #
#---------------------------------------------------------------------------------------#
znice = -pretty.log(-slz,n=8)
znice = sort(c(znice,slz[1],slz[nzg]))
sel = znice >= slz[1] & znice <= slz[nzg]
znice = znice[sel]
zat = -log(-znice)
nznice = length(znice)
znice = sprintf("%.2f",znice)
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#
# Define a suitable scale for monthly means... #
#---------------------------------------------------------------------------------------#
montmont = seq(from=1,to=12,by=1)
mplot = list()
mplot$levels = montmont
mplot$labels = capwords(mon2mmm(montmont))
mplot$n = 12
mplot$scale = "months"
mplot$padj = rep(0,times=mplot$n)
#---------------------------------------------------------------------------------------#
#=======================================================================================#
#=======================================================================================#
#=======================================================================================#
# Plotting section begins here... #
#---------------------------------------------------------------------------------------#
cat (" - Plotting figures...","\n")
#---------------------------------------------------------------------------------------#
#
# #---------------------------------------------------------------------------------------#
# # Time series by PFT. (tspft folder) #
# #---------------------------------------------------------------------------------------#
# for (v in sequence(ntspftdbh)){
# thistspft = tspftdbh[[v]]
# vnam = thistspft$vnam
# description = thistspft$desc
# unit = thistspft$e.unit
# plog = thistspft$plog
# plotit = thistspft$pft
#
# #----- Check whether the user wants to have this variable plotted. ------------------#
# if (plotit && any(selpft)){
#
# #---------------------------------------------------------------------------------#
# # Check whether the time series directory exists. If not, create it. #
# #---------------------------------------------------------------------------------#
# outdir = file.path(outpref,"tspft")
# if (! file.exists(outdir)) dir.create(outdir)
# cat(" +",description,"time series for all PFTs...","\n")
#
# #----- Load variable -------------------------------------------------------------#
# if (vnam %in% names(szpft)){
# thisvar = szpft[[vnam]][,ndbh+1,]
# if (plog){
# #----- Eliminate non-positive values in case it is a log plot. -------------#
# badlog = is.finite(thisvar) & thisvar <= 0
# thisvar[badlog] = NA
# }#end if
# }else{
# thisvar = matrix(NA,ncol=npft+1,nrow=nyears)
# }#end if
# #---------------------------------------------------------------------------------#
#
#
#
# #----- Loop over output formats. -------------------------------------------------#
# for (o in sequence(nout)){
# fichier = file.path(outdir,paste0(vnam,"-",suffix,".",outform[o]))
# if(outform[o] %in% "x11"){
# X11(width=size$width,height=size$height,pointsize=ptsz)
# }else if(outform[o] %in% "quartz"){
# quartz(width=size$width,height=size$height,pointsize=ptsz)
# }else if(outform[o] %in% "png"){
# png(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent")
# }else if(outform[o] %in% "tif"){
# tiff(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent",compression="lzw")
# }else if(outform[o] %in% "eps"){
# postscript(file=fichier,width=size$width,height=size$height
# ,pointsize=ptsz,paper=size$paper)
# }else if(outform[o] %in% "pdf"){
# pdf(file=fichier,onefile=FALSE
# ,width=size$width,height=size$height,pointsize=ptsz,paper=size$paper)
# }#end if
#
#
# #------------------------------------------------------------------------------#
# # Find the limit, make some room for the legend, and in case the field is #
# # a constant, nudge the limits so the plot command will not complain. #
# #------------------------------------------------------------------------------#
# xlimit = pretty.xylim(u = datum$toyear ,fracexp=0.0,is.log=FALSE)
# #manfredo 574 instead of ylimit = pretty.xylim(u = thisvar[,selpft],fracexp=0.0,is.log=plog )
# ylimit = pretty.xylim(u = thisvar[selpft],fracexp=0.2,is.log=plog )
# if (plog){
# xylog = "y"
# ydrought = c( exp(sqrt(ylimit[1]^3/ylimit[2]))
# , exp(sqrt(ylimit[2]^3/ylimit[1]))
# )#end c
# }else{
# xylog = ""
# ydrought = c( ylimit[1] - 0.5 * diff(ylimit),ylimit[2] + 0.5 * diff(ylimit) )
# }#end if
# #------------------------------------------------------------------------------#
#
#
# #----- Plot settings. ---------------------------------------------------------#
# letitre = paste(description,lieu,sep=" - ")
# ley = desc.unit(desc=description,unit=unit)
# cols = pft$colour[selpft]
# legs = pft$name [selpft]
# #------------------------------------------------------------------------------#
#
#
# #------------------------------------------------------------------------------#
# # Split the plot into two windows. #
# #------------------------------------------------------------------------------#
# par(par.user)
# layout(mat=rbind(2,1),heights=c(5,1))
# #------------------------------------------------------------------------------#
#
#
#
# #------------------------------------------------------------------------------#
# # First plot: legend. #
# #------------------------------------------------------------------------------#
# par(mar=c(0.1,4.6,0.1,2.1))
# plot.new()
# plot.window(xlim=c(0,1),ylim=c(0,1))
# legend( x = "bottom"
# , inset = 0.0
# , legend = legs
# , col = cols
# , lwd = lwidth
# , ncol = min(pretty.box(n.selpft)$ncol,3)
# , title = expression(bold("Plant Functional Type"))
# , xpd = TRUE
# , bty = "n"
# )#end legend
# #------------------------------------------------------------------------------#
#
#
#
# #------------------------------------------------------------------------------#
# # Main plot. #
# #------------------------------------------------------------------------------#
# par(mar=c(4.1,4.6,4.1,2.1))
# plot.new()
# plot.window(xlim=xlimit,ylim=ylimit,log=xylog)
# axis(side=1)
# axis(side=2,las=1)
# box()
# #title(main=letitre,xlab="Year",ylab=ley,cex.main=0.7,log=xylog)
# title(main=letitre,xlab="Year",ylab=ley,cex.main=0.7)
# if (drought.mark){
# for (n in sequence(ndrought)){
# rect(xleft = drought[[n]][1],ybottom = ydrought[1]
# ,xright = drought[[n]][2],ytop = ydrought[2]
# ,col = grid.colour,border=NA)
# }#end for
# }#end if
# #----- Plot grid. -------------------------------------------------------------#
# if (plotgrid){
# abline(v=axTicks(side=1),h=axTicks(side=2),col=grid.colour,lty="solid")
# }#end if
# #----- Plot lines. ------------------------------------------------------------#
# for (n in sequence(npft+1)){
# if (selpft[n]){
# lines(datum$toyear,thisvar[,n],type="l",col=pft$colour[n],lwd=lwidth)
# }#end if
# }#end for
# #------------------------------------------------------------------------------#
#
#
# #----- Close the device. ------------------------------------------------------#
# if (outform[o] %in% c("x11","quartz")){
# locator(n=1)
# dev.off()
# }else{
# dev.off()
# }#end if
# dummy=clean.tmp()
# #------------------------------------------------------------------------------#
# } #end for outform
# }#end if (tseragbpft)
# } #end for tseries
# #---------------------------------------------------------------------------------------#
#
#
#
#
# #---------------------------------------------------------------------------------------#
# # Time series by DBH, by PFT. (tsdbh folder) #
# #---------------------------------------------------------------------------------------#
# #----- Find the PFTs to plot. ----------------------------------------------------------#
# pftuse = which(apply(X=is.na(szpft$nplant),MARGIN=3,FUN=sum,na.rm=TRUE) == 0.)
# pftuse = pftuse[pftuse != (npft+1)]
# for (v in sequence(ntspftdbh)){
# thistspftdbh = tspftdbh[[v]]
# vnam = thistspftdbh$vnam
# description = thistspftdbh$desc
# unit = thistspftdbh$e.unit
# plog = thistspftdbh$plog
# plotit = thistspftdbh$pftdbh
#
# #----- Load variable ----------------------------------------------------------------#
# if (vnam %in% names(szpft)){
# thisvar = szpft[[vnam]]
# if (plog){
# xylog="y"
# badlog = is.finite(thisvar) & thisvar <= 0
# thisvar[badlog] = NA
# }else{
# xylog=""
# }#end if
# }else{
# thisvar = array(NA,dim=c(nyears,ndbh+1,npft+1))
# }#end if
# #----- Check whether the user wants to have this variable plotted. ------------------#
# if (plotit && length(pftuse) > 0 && any(is.finite(thisvar))){
#
# #---------------------------------------------------------------------------------#
# # Check whether the time series directory exists. If not, create it. #
# #---------------------------------------------------------------------------------#
# outdir = file.path(outpref,"tsdbh")
# if (! file.exists(outdir)) dir.create(outdir)
# outvar = file.path(outdir,vnam)
# if (! file.exists(outvar)) dir.create(outvar)
# #---------------------------------------------------------------------------------#
#
# cat(" +",description,"time series for DBH class...","\n")
#
#
# #---------------------------------------------------------------------------------#
# # Find the limit, make some room for the legend, and in case the field is a #
# # constant, nudge the limits so the plot command will not complain. #
# #---------------------------------------------------------------------------------#
# xlimit = pretty.xylim(u=datum$toyear ,fracexp=0.0,is.log=FALSE)
# ylimit = pretty.xylim(u=thisvar[,,pftuse],fracexp=0.0,is.log=plog)
# if (plog){
# xylog = "y"
# ydrought = c( exp(sqrt(ylimit[1]^3/ylimit[2]))
# , exp(sqrt(ylimit[2]^3/ylimit[1]))
# )#end c
# }else{
# xylog = ""
# ydrought = c( ylimit[1] - 0.5 * diff(ylimit),ylimit[2] + 0.5 * diff(ylimit) )
# }#end if
# #---------------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------------#
# # Loop over plant functional types. #
# #---------------------------------------------------------------------------------#
# for (p in pftuse){
# pftlab = paste0("pft-",sprintf("%2.2i",p))
# cat(" - ",pft$name[p],"\n")
#
#
# #----- Loop over output formats. ----------------------------------------------#
# for (o in sequence(nout)){
# #----- Open file. ----------------------------------------------------------#
# fichier = file.path( outvar
# , paste0(vnam,"-",pftlab,"-",suffix,".",outform[o])
# )#end file.path
# if (outform[o] %in% "x11"){
# X11(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "quartz"){
# quartz(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "png"){
# png(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent")
# }else if (outform[o] %in% "tif"){
# tiff(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent",compression="lzw")
# }else if (outform[o] %in% "eps"){
# postscript(file=fichier,width=size$width,height=size$height
# ,pointsize=ptsz,paper=size$paper)
# }else if (outform[o] %in% "pdf"){
# pdf(file=fichier,onefile=FALSE
# ,width=size$width,height=size$height,pointsize=ptsz,paper=size$paper)
# }#end if
# #---------------------------------------------------------------------------#
#
#
#
# #----- Plot annotation. ---------------------------------------------------#
# letitre = paste(description,pft$name[p],lieu,sep=" - ")
# ley = desc.unit(desc=description,unit=unit)
# #---------------------------------------------------------------------------#
#
#
# #---------------------------------------------------------------------------#
# # Split the plot into two windows. #
# #---------------------------------------------------------------------------#
# par(par.user)
# layout(mat=rbind(2,1),heights=c(5,1))
# #---------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------#
# # First plot: legend. #
# #---------------------------------------------------------------------------#
# par(mar=c(0.1,4.6,0.1,2.1))
# plot.new()
# plot.window(xlim=c(0,1),ylim=c(0,1))
# legend( x = "bottom"
# , inset = 0.0
# , bg = background
# , legend = dbhnames
# , col = dbhcols
# , ncol = min(pretty.box(ndbh+1)$ncol,3)
# , title = expression(bold("DBH class"))
# , lwd = lwidth
# , xpd = TRUE
# , bty = "n"
# )#end legend
# #---------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------#
# # Main plot. #
# #---------------------------------------------------------------------------#
# par(mar=c(4.1,4.6,4.1,2.1))
# plot.new()
# plot.window(xlim=xlimit,ylim=ylimit,log=xylog)
# axis(side=1)
# axis(side=2,las=1)
# box()
# title(main=letitre,xlab="Year",ylab=ley,cex.main=0.7)
# if (drought.mark){
# for (n in sequence(ndrought)){
# rect(xleft = drought[[n]][1],ybottom = ydrought[1]
# ,xright = drought[[n]][2],ytop = ydrought[2]
# ,col = grid.colour,border=NA)
# }#end for
# }#end if
# #----- Plot grid. ----------------------------------------------------------#
# if (plotgrid){
# abline(v=axTicks(side=1),h=axTicks(side=2),col=grid.colour,lty="solid")
# }#end if
# #----- Plot lines. ---------------------------------------------------------#
# for (d in seq(from=1,to=ndbh+1,by=1)){
# lines(datum$toyear,thisvar[,d,p],type="l",col=dbhcols[d],lwd=lwidth)
# }#end for
# #---------------------------------------------------------------------------#
#
#
#
# #----- Close the device. ---------------------------------------------------#
# if (outform[o] %in% c("x11","quartz")){
# locator(n=1)
# dev.off()
# }else{
# dev.off()
# }#end if
# dummy=clean.tmp()
# #---------------------------------------------------------------------------#
# }#end for outform
# #------------------------------------------------------------------------------#
# }#end for (p in pftuse)
# #---------------------------------------------------------------------------------#
# }#end if (tseragbpft)
# #------------------------------------------------------------------------------------#
# } #end for tseries
# #---------------------------------------------------------------------------------------#
#
#
#
#
# #---------------------------------------------------------------------------------------#
# # Plot the comparison between observations and model. (ycomp folder) #
# #---------------------------------------------------------------------------------------#
# cat(" + Year-by-year comparisons of monthly means...","\n")
# for (cc in sequence(ncompmodel)){
#
# #----- Retrieve variable information from the list. ---------------------------------#
# compnow = compmodel[[cc]]
# vname = compnow$vnam
# description = compnow$desc
# unit = compnow$unit
# plotsd = compnow$plotsd
# lcolours = compnow$colour
# errcolours = compnow$errcol
# angle = compnow$angle
# dens = compnow$dens
# llwd = compnow$lwd
# shwd = compnow$shwd
# ltype = compnow$type
# plog = compnow$plog
# legpos = compnow$legpos
# plotit = compnow$mmean
#
# plotit = ( plotit && vname %in% names(emean) && vname %in% names(mmean)
# && plot.ycomp )
#
# if (plotit){
# #---------------------------------------------------------------------------------#
# # Copy the observations to a scratch variable. #
# #---------------------------------------------------------------------------------#
# thisvar = emean [[vname]]
# thismean = mmean [[vname]]
# if (length(msdev[[vname]]) == 0){
# thissdev = 0. * thismean
# }else{
# thissdev = msdev[[vname]]
# }#end if
# mod.x = montmont
# mod.ylow = thismean - thissdev
# mod.yhigh = thismean + thissdev
# mod.x.poly = c(mod.x,rev(mod.x))
# mod.y.poly = c(mod.ylow,rev(mod.yhigh))
# mod.keep = is.finite(mod.y.poly)
# mod.x.poly = mod.x.poly[mod.keep]
# mod.y.poly = mod.y.poly[mod.keep]
# #---------------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------------#
# # Check whether the time series directory exists. If not, create it. #
# #---------------------------------------------------------------------------------#
# outdir = file.path(outpref,"ycomp")
# outvar = file.path(outdir,vname)
# if (! file.exists(outdir)) dir.create(outdir)
# if (! file.exists(outvar)) dir.create(outvar)
# cat(" - ",description,"comparison...","\n")
# #---------------------------------------------------------------------------------#
#
#
#
# #----- Find the plot range. ------------------------------------------------------#
# if (plotsd){
# ylimit = range(c(mod.ylow,mod.yhigh,thisvar),na.rm=TRUE)
# }else{
# ylimit = range(thisvar,na.rm=TRUE)
# }#end if
# #----- Expand the upper range in so the legend doesn't hide things. --------------#
# ylimit = pretty.xylim(u=ylimit,fracexp=0.0,is.log=FALSE)
# #---------------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------------#
# # Loop over all years, and make one plot per year. #
# #---------------------------------------------------------------------------------#
# for (y in sequence(nyears)){
# #----- Retrieve the year and the variable for this year. ----------------------#
# year.now = datum$toyear[y]
# cyear = sprintf("%4.4i",year.now)
# var.year = thisvar[yfac == year.now]
# #------------------------------------------------------------------------------#
#
#
#
# #----- Load variable ----------------------------------------------------------#
# letitre = paste0(description," - ",lieu,"\n","Monthly mean - ",cyear)
# ley = desc.unit(desc=description,unit=unit)
# #------------------------------------------------------------------------------#
#
#
# #----- Loop over formats. -----------------------------------------------------#
# for (o in sequence(nout)){
# fichier = file.path(outvar,paste0(vname,"-",cyear,".",outform[o]))
# if (outform[o] %in% "x11"){
# X11(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "quartz"){
# quartz(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "png"){
# png(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent")
# }else if (outform[o] %in% "tif"){
# tiff(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent",compression="lzw")
# }else if (outform[o] %in% "eps"){
# postscript(file=fichier,width=size$width,height=size$height
# ,pointsize=ptsz,paper=size$paper)
# }else if (outform[o] %in% "pdf"){
# pdf(file=fichier,onefile=FALSE
# ,width=size$width,height=size$height,pointsize=ptsz,paper=size$paper)
# }#end if
#
#
#
# #----- Split plot into two windows. ----------------------------------------#
# par(par.user)
# layout(mat=rbind(2,1),heights=c(5,1))
# #---------------------------------------------------------------------------#
#
#
#
# #------ First plot: the legend. --------------------------------------------#
# par(mar=c(0.1,4.6,0.1,2.1))
# plot.new()
# plot.window(xlim=c(0,1),ylim=c(0,1))
# if (plotsd){
# legend( x = "bottom"
# , inset = 0.0
# , legend = c(cyear,paste0("Mean: ",yeara,"-",yearz))
# , fill = errcolours
# , angle = angle
# , density = dens
# , lwd = llwd
# , col = lcolours
# , bg = background
# , title = expression(bold("Shaded areas = 1 SD"))
# , cex = cex.ptsz
# , pch = 16
# , xpd = TRUE
# , bty = "n"
# )#end legend
# }else{
# legend( x = "bottom"
# , inset = 0.0
# , legend = c(cyear,paste0("Mean: ",yeara,"-",yearz))
# , col = lcolours
# , lwd = llwd
# , cex = cex.ptsz
# , pch = 16
# , xpd = TRUE
# , bty = "n"
# )#end legend
# }#end if
# #---------------------------------------------------------------------------#
#
#
#
#
# #------ Second plot: the comparison. ---------------------------------------#
# par(mar=c(4.1,4.6,4.1,2.1))
# plot.new()
# plot.window(xlim=range(montmont),ylim=ylimit,log=plog)
# if (plotgrid){
# abline(v=mplot$levels,h=axTicks(side=2),col=grid.colour,lty="solid")
# }#end if
# if (plotsd){
# polygon(x=mod.x.poly,y=mod.y.poly,col=errcolours[2],angle=angle[2]
# ,density=dens[1],lty="solid",lwd=shwd[1])
# }#end if
# points(x=montmont,y=var.year,col=lcolours[1],lwd=llwd[1],type=ltype
# ,pch=16,cex=1.0)
# points(x=montmont,y=thismean,col=lcolours[2],lwd=llwd[2],type=ltype
# ,pch=16,cex=1.0)
# axis(side=1,at=mplot$levels,labels=mplot$labels,padj=mplot$padj)
# axis(side=2,las=1)
# title(main=letitre,xlab="Time",ylab=ley,cex.main=cex.main)
# box()
# #---------------------------------------------------------------------------#
#
#
#
#
# #----- Close plotting window. ----------------------------------------------#
# if (outform[o] %in% c("x11","quartz")){
# locator(n=1)
# dev.off()
# }else{
# dev.off()
# }#end if
# dummy=clean.tmp()
# #---------------------------------------------------------------------------#
# }#end for outform
# #------------------------------------------------------------------------------#
# }#end for years
# #---------------------------------------------------------------------------------#
# }#end if plotit
# #------------------------------------------------------------------------------------#
# }#end for ncompare
# #---------------------------------------------------------------------------------------#
#
#
#
#
# #---------------------------------------------------------------------------------------#
# # Time series by LU. (tslu folder) #
# #---------------------------------------------------------------------------------------#
# for (v in sequence(ntslu)){
# thistslu = tslu[[v]]
# vnam = thistslu$vnam
# description = thistslu$desc
# unit = thistslu$unit
# plog = thistslu$plog
# plotit = thistslu$plt
#
# #----- Check whether the user wants to have this variable plotted. ------------------#
# if (plotit && any(sellu)){
#
# #---------------------------------------------------------------------------------#
# # Check whether the time series directory exists. If not, create it. #
# #---------------------------------------------------------------------------------#
# outdir = file.path(outpref,"tslu")
# if (! file.exists(outdir)) dir.create(outdir)
# cat(" +",description,"time series for all LUs...","\n")
#
#
#
# #----- Load variable -------------------------------------------------------------#
# thisvar = lu[[vnam]]
# if (plog){
# #----- Eliminate non-positive values in case it is a log plot. ----------------#
# thisvar[thisvar <= 0] = NA
# }#end if
# #---------------------------------------------------------------------------------#
#
# #----- Loop over output formats. -------------------------------------------------#
# for (o in sequence(nout)){
# fichier = file.path(outdir,paste0(vnam,"-",suffix,".",outform[o]))
# if (outform[o] %in% "x11"){
# X11(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "quartz"){
# quartz(width=size$width,height=size$height,pointsize=ptsz)
# }else if(outform[o] %in% "png"){
# png(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent")
# }else if(outform[o] %in% "tif"){
# tiff(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent",compression="lzw")
# }else if(outform[o] %in% "eps"){
# postscript(file=fichier,width=size$width,height=size$height
# ,pointsize=ptsz,paper=size$paper)
# }else if(outform[o] %in% "pdf"){
# pdf(file=fichier,onefile=FALSE
# ,width=size$width,height=size$height,pointsize=ptsz,paper=size$paper)
# }#end if
#
#
# #------------------------------------------------------------------------------#
# # Find the limit, make some room for the legend, and in case the field is #
# # a constant, nudge the limits so the plot command will not complain. #
# #------------------------------------------------------------------------------#
# xlimit = pretty.xylim(u = datum$toyear ,fracexp=0.0,is.log=FALSE)
# ylimit = pretty.xylim(u=thisvar[,sellu],fracexp=0.0,is.log=plog)
# if (plog){
# xylog = "y"
# ydrought = c( exp(sqrt(ylimit[1]^3/ylimit[2]))
# , exp(sqrt(ylimit[2]^3/ylimit[1]))
# )#end c
# }else{
# xylog = ""
# ydrought = c( ylimit[1] - 0.5 * diff(ylimit),ylimit[2] + 0.5 * diff(ylimit) )
# }#end if
# #------------------------------------------------------------------------------#
#
#
#
# #----- Plot settings. ---------------------------------------------------------#
# letitre = paste(description,lieu,sep=" - ")
# ley = desc.unit(desc=description,unit=unit)
# cols = lucols[sellu]
# legs = lunames[sellu]
# #------------------------------------------------------------------------------#
#
#
# #------------------------------------------------------------------------------#
# # Split the plot into two windows. #
# #------------------------------------------------------------------------------#
# par(par.user)
# layout(mat=rbind(2,1),heights=c(5,1))
# #------------------------------------------------------------------------------#
#
#
#
# #------------------------------------------------------------------------------#
# # First plot: legend. #
# #------------------------------------------------------------------------------#
# par(mar=c(0.1,4.6,0.1,2.1))
# plot.new()
# plot.window(xlim=c(0,1),ylim=c(0,1))
# legend( x = "bottom"
# , inset = 0.0
# , legend = legs
# , col = cols
# , lwd = lwidth
# , ncol = min(3,pretty.box(n.sellu)$ncol)
# , title = expression(bold("Land use type"))
# , xpd = TRUE
# , bty = "n"
# )#end legend
# #------------------------------------------------------------------------------#
#
#
#
# #------------------------------------------------------------------------------#
# # Main plot. #
# #------------------------------------------------------------------------------#
# par(mar=c(4.1,4.6,4.1,2.1))
# plot.new()
# plot.window(xlim=xlimit,ylim=ylimit,log=xylog)
# axis(side=1)
# axis(side=2,las=1)
# box()
# title(main=letitre,xlab="Year",ylab=ley,cex.main=0.7)
# if (drought.mark){
# for (n in sequence(ndrought)){
# rect(xleft = drought[[n]][1],ybottom = ydrought[1]
# ,xright = drought[[n]][2],ytop = ydrought[2]
# ,col = grid.colour,border=NA)
# }#end for
# }#end if
# #----- Plot grid. -------------------------------------------------------------#
# if (plotgrid){
# abline(v=axTicks(side=1),h=axTicks(side=2),col=grid.colour,lty="solid")
# }#end if
# #----- Plot lines. ------------------------------------------------------------#
# for (n in sequence(nlu+1)){
# if (sellu[n]){
# lines(datum$toyear,thisvar[,n],type="l",col=lucols[n],lwd=lwidth)
# }#end if
# }#end for
# #------------------------------------------------------------------------------#
#
#
# #----- Close the device. ------------------------------------------------------#
# if (outform[o] %in% c("x11","quartz")){
# locator(n=1)
# dev.off()
# }else{
# dev.off()
# }#end if
# dummy=clean.tmp()
# #------------------------------------------------------------------------------#
# }#end for outform
# #---------------------------------------------------------------------------------#
# }#end if (tseragbpft)
# #------------------------------------------------------------------------------------#
# }#end for tseries
# #---------------------------------------------------------------------------------------#
#
#
#
#
#
# #---------------------------------------------------------------------------------------#
# # Plot disturbance rate by disturbance transition. (distrurb-gyf.pdf file) #
# #---------------------------------------------------------------------------------------#
# if (tserdist && any(seldist)){
# cat(" + Disturbance rate time series for all disturbances...","\n")
# for (o in sequence(nout)){
# fichier = file.path(outpref,paste0("disturb-",suffix,".",outform[o]))
# if (outform[o] %in% "x11"){
# X11(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "quartz"){
# quartz(width=size$width,height=size$height,pointsize=ptsz)
# }else if(outform[o] %in% "png"){
# png(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent")
# }else if(outform[o] %in% "tif"){
# tiff(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent",compression="lzw")
# }else if(outform[o] %in% "eps"){
# postscript(file=fichier,width=size$width,height=size$height
# ,pointsize=ptsz,paper=size$paper)
# }else if(outform[o] %in% "pdf"){
# pdf(file=fichier,onefile=FALSE
# ,width=size$width,height=size$height,pointsize=ptsz,paper=size$paper)
# }#end if
#
# #---------------------------------------------------------------------------------#
# # Find the limit, make some room for the legend, and in case the field is a #
# # constant, nudge the limits so the plot command will not complain. #
# #---------------------------------------------------------------------------------#
# xlimit = pretty.xylim(u=datum$toyear,fracexp=0.0,is.log=FALSE)
# ylimit = NULL
# n = 0
# mylucols = NULL
# mylulegs = NULL
# for (jlu in sequence(nlu)){
# for (ilu in sequence(nlu)){
# n = n + 1
# if (seldist[ilu,jlu]){
# ylimit = c(ylimit,lu$dist[,ilu,jlu])
# mylucols = c(mylucols,distcols [n])
# mylulegs = c(mylulegs,distnames[n])
# }#end if
# }#end for
# }#end for
# ylimit = pretty.xylim(u=ylimit,fracexp=0.0,is.log=FALSE)
# ydrought = c(ylimit[1] - 0.5 * diff(ylimit), ylimit[2] + 0.5 * diff(ylimit))
# #---------------------------------------------------------------------------------#
#
#
#
# #----- Plot settings. ------------------------------------------------------------#
# letitre = paste("Disturbance rates",lieu,sep=" - ")
# #---------------------------------------------------------------------------------#
#
#
# #---------------------------------------------------------------------------------#
# # Split the plot into two windows. #
# #---------------------------------------------------------------------------------#
# par(par.user)
# layout(mat=rbind(2,1),heights=c(5,1))
# #---------------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------------#
# # First plot: legend. #
# #---------------------------------------------------------------------------------#
# par(mar=c(0.1,4.6,0.1,2.1))
# plot.new()
# plot.window(xlim=c(0,1),ylim=c(0,1))
# legend( x = "bottom"
# , inset = 0.0
# , bg = background
# , legend = mylulegs
# , col = mylucols
# , lwd = lwidth
# , ncol = min(3,pretty.box(n)$ncol)
# , title = expression(bold("Transition"))
# , xpd = TRUE
# , bty = "n"
# )#end legend
# #---------------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------------#
# # Main plot. #
# #---------------------------------------------------------------------------------#
# par(mar=c(4.1,4.6,4.1,2.1))
# plot.new()
# plot.window(xlim=xlimit,ylim=ylimit,log=xylog)
# axis(side=1)
# axis(side=2,las=1)
# box()
# title( main = letitre
# , xlab = "Year"
# , ylab = desc.unit(desc="Disturbance rate",unit=untab$oneoyr)
# , cex.main = 0.7
# )#end title
# if (drought.mark){
# for (n in sequence(ndrought)){
# rect(xleft = drought[[n]][1],ybottom = ydrought[1]
# ,xright = drought[[n]][2],ytop = ydrought[2]
# ,col = grid.colour,border=NA)
# }#end for
# }#end if
# #----- Plot grid. ----------------------------------------------------------------#
# if (plotgrid){
# abline(v=axTicks(side=1),h=axTicks(side=2),col=grid.colour,lty="solid")
# }#end if
# #----- Plot lines. ---------------------------------------------------------------#
# n = 0
# for (jlu in sequence(nlu)){
# for (ilu in sequence(nlu)){
# n = n + 1
# if (seldist[ilu,jlu]){
# lines(datum$toyear,lu$dist[,ilu,jlu],type="l"
# ,col=distcols[n],lwd=lwidth)
# }#end if
# }#end for
# }#end for
# #---------------------------------------------------------------------------------#
#
#
# #----- Close the device. ---------------------------------------------------------#
# if (outform[o] %in% c("x11","quartz")){
# locator(n=1)
# dev.off()
# }else{
# dev.off()
# }#end if
# dummy=clean.tmp()
# #---------------------------------------------------------------------------------#
# } #end for outform
# #------------------------------------------------------------------------------------#
# }#end if
# #---------------------------------------------------------------------------------------#
#
#
#
#
# #---------------------------------------------------------------------------------------#
# # Plot the time series diagrams showing annual means. (theme_ymean folder) #
# #---------------------------------------------------------------------------------------#
# cat(" * Plot time series of groups of variables...","\n")
# for (hh in sequence(ntheme)){
#
# #----- Retrieve variable information from the list. ---------------------------------#
# themenow = theme[[hh]]
# vnames = themenow$vnam
# description = themenow$desc
# lcolours = themenow$colour
# llwd = themenow$lwd
# ltype = themenow$type
# plog = themenow$plog
# prefix = themenow$prefix
# group = themenow$title
# unit = themenow$unit
# legpos = themenow$legpos
# plotit = themenow$ymean
# ylimit.fix = themenow$ymean.lim
#
# if (plotit){
#
# #---------------------------------------------------------------------------------#
# # Check whether the time series directory exists. If not, create it. #
# #---------------------------------------------------------------------------------#
# outdir = file.path(outpref,"theme_ymean")
# if (! file.exists(outdir)) dir.create(outdir)
# cat(" +",group,"time series...","\n")
#
#
# #----- Define the number of layers. ----------------------------------------------#
# nlayers = length(vnames)
# #---------------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------------#
# # Find the limit, make some room for the legend, and in case the field is a #
# # constant, nudge the limits so the plot command will not complain. #
# #---------------------------------------------------------------------------------#
# xlimit = pretty.xylim(u=datum$toyear,fracexp=0.0,is.log=FALSE)
# if (any(! is.finite(ylimit.fix))){
# ylimit = NULL
# for (l in sequence(nlayers)){
# thisvar = ymean[[vnames[l]]]
# ylimit = range(c(ylimit,thisvar),na.rm=TRUE)
# }#end for
# ylimit = pretty.xylim(u=ylimit,fracexp=0.0,is.log=plog)
# }else{
# ylimit = ylimit.fix
# }#end if
# if (plog) {
# xylog = "y"
# ydrought = c( exp(sqrt(ylimit[1]^3/ylimit[2]))
# , exp(sqrt(ylimit[2]^3/ylimit[1]))
# )#end c
# }else{
# xylog = ""
# ydrought = c(ylimit[1]-0.5*diff(ylimit),ylimit[2]+0.5*diff(ylimit))
# }#end if
# #---------------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------------#
# # Check if the directory exists. If not, create it. #
# #---------------------------------------------------------------------------------#
#
# #----- Loop over formats. --------------------------------------------------------#
# for (o in sequence(nout)){
# #------ Open file. ------------------------------------------------------------#
# fichier = file.path(outdir,paste0(prefix,"-",suffix,".",outform[o]))
# if (outform[o] %in% "x11"){
# X11(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "quartz"){
# quartz(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "png"){
# png(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent")
# }else if (outform[o] %in% "tiff"){
# tiff(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent",compression="lzw")
# }else if (outform[o] %in% "eps"){
# postscript(file=fichier,width=size$width,height=size$height
# ,pointsize=ptsz,paper=size$paper)
# }else if (outform[o] %in% "pdf"){
# pdf(file=fichier,onefile=FALSE
# ,width=size$width,height=size$height,pointsize=ptsz,paper=size$paper)
# }#end if
# #------------------------------------------------------------------------------#
#
#
#
# #----- Plot settings. ---------------------------------------------------------#
# letitre = paste0(" Time series: ",group,"\n",lieu)
# ley = desc.unit(desc=group,unit=unit)
# #------------------------------------------------------------------------------#
#
#
# #------------------------------------------------------------------------------#
# # Split the plot into two windows. #
# #------------------------------------------------------------------------------#
# par(par.user)
# layout(mat=rbind(2,1),heights=c(5,1))
# #------------------------------------------------------------------------------#
#
#
#
# #------------------------------------------------------------------------------#
# # First plot: legend. #
# #------------------------------------------------------------------------------#
# par(mar=c(0.1,4.6,0.1,2.1))
# plot.new()
# plot.window(xlim=c(0,1),ylim=c(0,1))
# legend( x = "bottom"
# , inset = 0.0
# , legend = description
# , col = lcolours
# , lwd = llwd
# , ncol = min(3,pretty.box(nlayers)$ncol)
# , xpd = TRUE
# , bty = "n"
# )#end legend
# #------------------------------------------------------------------------------#
#
#
#
# #------------------------------------------------------------------------------#
# # Main plot. #
# #------------------------------------------------------------------------------#
# par(mar=c(4.1,4.6,4.1,2.1))
# plot.new()
# plot.window(xlim=xlimit,ylim=ylimit,log=xylog)
# axis(side=1)
# axis(side=2,las=1)
# box()
# title(main=letitre,xlab="Year",ylab=ley,cex.main=0.7)
# if (drought.mark){
# for (n in sequence(ndrought)){
# rect(xleft = drought[[n]][1],ybottom = ydrought[1]
# ,xright = drought[[n]][2],ytop = ydrought[2]
# ,col = grid.colour,border=NA)
# }#end for
# }#end if
# #----- Plot grid. -------------------------------------------------------------#
# if (plotgrid){
# abline(v=axTicks(side=1),h=axTicks(side=2),col=grid.colour,lty="solid")
# }#end if
# #----- Plot lines. ------------------------------------------------------------#
# for (l in sequence(nlayers)){
# thisvar = ymean[[vnames[l]]]
# points(x=datum$toyear,y=thisvar,col=lcolours[l],lwd=llwd[l],type=ltype
# ,pch=16,cex=0.8)
# }#end for
# #------------------------------------------------------------------------------#
#
#
# #----- Close the device. ------------------------------------------------------#
# if (outform[o] %in% c("x11","quartz")){
# locator(n=1)
# dev.off()
# }else{
# dev.off()
# }#end if
# dummy=clean.tmp()
# #------------------------------------------------------------------------------#
# } #end for outform
# #---------------------------------------------------------------------------------#
# }#end if plotit
# #------------------------------------------------------------------------------------#
# }#end for ntser
# #---------------------------------------------------------------------------------------#
#
#
#
#
# #---------------------------------------------------------------------------------------#
# # Plot the climatology of the mean diurnal cycle. (theme_qmean folder) #
# #---------------------------------------------------------------------------------------#
# cat(" * Plot mean diel for groups of variables...","\n")
# for (hh in sequence(ntheme)){
#
# #----- Retrieve variable information from the list. ---------------------------------#
# themenow = theme[[hh]]
# vnames = themenow$vnam
# description = themenow$desc
# lcolours = themenow$colour
# llwd = themenow$lwd
# ltype = themenow$type
# plog = themenow$plog
# prefix = themenow$prefix
# group = themenow$title
# unit = themenow$unit
# legpos = themenow$legpos
# ylimit.fix = themenow$qmean.lim
# plotit = themenow$qmean && plot.ycomp
# if (plog){
# xylog = "y"
# }else{
# xylog = ""
# }#end if
#
# if (plotit){
#
# #---------------------------------------------------------------------------------#
# # Check whether the time series directory exists. If not, create it. #
# #---------------------------------------------------------------------------------#
# outdir = file.path(outpref,"theme_qmean")
# if (! file.exists(outdir)) dir.create(outdir)
# outtheme = file.path(outdir,prefix)
# if (! file.exists(outtheme)) dir.create(outtheme)
# cat(" +",group," diurnal cycle...","\n")
#
#
# #----- Define the number of layers. ----------------------------------------------#
# nlayers = length(vnames)
# xlimit = range(thisday)
# if (any(! is.finite(ylimit.fix))){
# ylimit = NULL
# for (l in sequence(nlayers)) ylimit = c(ylimit,umean[[vnames[l]]])
# ylimit = pretty.xylim(u=ylimit,fracexp=0.0,is.log=plog)
# }else{
# ylimit = ylimit.fix
# }#end if
# #---------------------------------------------------------------------------------#
#
#
# #---------------------------------------------------------------------------------#
# # Loop over all months. #
# #---------------------------------------------------------------------------------#
# yplot = as.numeric(dimnames(umean[[vnames[1]]])[[1]])
# for (yy in seq_along(yplot)){
# cyear = sprintf("%4.4i",yplot[yy])
#
# #------------------------------------------------------------------------------#
# # Check if the directory exists. If not, create it. #
# #------------------------------------------------------------------------------#
#
# #----- Loop over formats. -----------------------------------------------------#
# for (o in sequence(nout)){
# #------ Open file. ---------------------------------------------------------#
# fichier = file.path( outtheme
# , paste0(prefix,"-",cyear,"-",suffix,".",outform[o])
# )#end file.path
# if (outform[o] %in% "x11"){
# X11(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "quartz"){
# quartz(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "png"){
# png(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent")
# }else if (outform[o] %in% "tif"){
# tiff(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent",compression="lzw")
# }else if (outform[o] %in% "eps"){
# postscript(file=fichier,width=size$width,height=size$height
# ,pointsize=ptsz,paper=size$paper)
# }else if (outform[o] %in% "pdf"){
# pdf(file=fichier,onefile=FALSE
# ,width=size$width,height=size$height,pointsize=ptsz,paper=size$paper)
# }#end if
# #---------------------------------------------------------------------------#
#
#
#
# #----- Plot settings. ------------------------------------------------------#
# letitre = paste0(group," - ",lieu,"\n","Mean diurnal cycle - ",cyear)
# ley = desc.unit(desc=group,unit=unit)
# #---------------------------------------------------------------------------#
#
#
# #---------------------------------------------------------------------------#
# # Split the plot into two windows. #
# #---------------------------------------------------------------------------#
# par(par.user)
# layout(mat=rbind(2,1),heights=c(5,1))
# #------------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------#
# # First plot: legend. #
# #---------------------------------------------------------------------------#
# par(mar=c(0.1,4.6,0.1,2.1))
# plot.new()
# plot.window(xlim=c(0,1),ylim=c(0,1))
# legend( x = "bottom"
# , inset = 0.0
# , legend = description
# , col = lcolours
# , lwd = llwd
# , ncol = min(3,pretty.box(nlayers)$ncol)
# , xpd = TRUE
# , bty = "n"
# )#end legend
# #---------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------#
# # Main plot. #
# #---------------------------------------------------------------------------#
# par(mar=c(4.1,4.6,4.1,2.1))
# plot.new()
# plot.window(xlim=xlimit,ylim=ylimit,log=xylog)
# axis(side=1,at=uplot$levels,labels=uplot$labels,padj=uplot$padj)
# axis(side=2,las=1)
# box()
# title(main=letitre,xlab="Year",ylab=ley,cex.main=0.7)
# #----- Plot grid. ----------------------------------------------------------#
# if (plotgrid){
# abline(v=uplot$levels,h=axTicks(side=2),col=grid.colour,lty="solid")
# }#end if
# #----- Plot lines. ---------------------------------------------------------#
# for (l in sequence(nlayers)){
# thisvar = umean[[vnames[l]]]
# thisvar = cbind(thisvar[,ndcycle],thisvar)
# points(x=thisday,y=thisvar[yy,],col=lcolours[l]
# ,lwd=llwd[l],type=ltype,pch=16)
# }#end for
# #---------------------------------------------------------------------------#
#
#
# #----- Close the device. ---------------------------------------------------#
# if (outform[o] %in% c("x11","quartz")){
# locator(n=1)
# dev.off()
# }else{
# dev.off()
# }#end if
# dummy=clean.tmp()
# #---------------------------------------------------------------------------#
# } #end for outform
# #------------------------------------------------------------------------------#
# }#end for pmon
# #---------------------------------------------------------------------------------#
# }#end if plotit
# #------------------------------------------------------------------------------------#
# }#end for ntser
# #---------------------------------------------------------------------------------------#
#
#
#
#
#
# #---------------------------------------------------------------------------------------#
# # Plot the climatology of the soil properties. (soil_ymean folder) #
# #---------------------------------------------------------------------------------------#
# for (v in sequence(nsoilplot)){
#
# #----- Retrieve variable information from the list. ---------------------------------#
# thisclim = soilplot[[v]]
# vnam = thisclim$vnam
# description = thisclim$desc
# unit = thisclim$unit
# vcscheme = thisclim$csch
# pnlog = thisclim$pnlog
# plotit = thisclim$ymean
#
# if (plotit){
#
# #---------------------------------------------------------------------------------#
# # Check if the directory exists. If not, create it. #
# #---------------------------------------------------------------------------------#
# outdir = file.path(outpref,"soil_ymean")
# if (! file.exists(outdir)) dir.create(outdir)
# cat(" + Climatology profile of ",description,"...","\n")
#
# #----- Find the number of rows and columns, and the axes. ------------------------#
# monaxis = sort(unique(datum$year))
# soilaxis = slz
# nmon = length(monaxis)
# nsoil = nzg
#
# #----- Convert the vector data into an array. ------------------------------------#
# vararr = ymean[[vnam]]
#
# #----- Copy the first and the last year to make the edges buffered. --------------#
# first = vararr[1,]
# first = c(first,first[nzg],first[nzg])
#
# last = vararr[nyears,]
# last = c(last[1],last[1],last)
# #---------------------------------------------------------------------------------#
#
#
#
# #----- Bind first and last year to the array, to make the edges buffered. --------#
# varbuff = cbind(vararr[,1],vararr,vararr[,nzg])
# varbuff = rbind(last,varbuff,first)
# #---------------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------------#
# # Expand the month and year axes. #
# #---------------------------------------------------------------------------------#
# yearaxis = c(min(datum$toyear)-1,datum$toyear,max(datum$toyear)+1)
# soilaxis = -log(-1.0 * c( slz[1]*(slz[1]/slz[2])
# , soilaxis
# , slz[nzg]*(slz[nzg]/slz[nzg-1]) ))
#
# if (pnlog){
# vrange = range(varbuff,na.rm=TRUE)
# vlevels = pretty.log(x=vrange,n=ncolshov)
# vnlev = length(vlevels)
# }else{
# vrange = range(varbuff,na.rm=TRUE)
# vlevels = pretty(x=vrange,n=ncolshov)
# vnlev = length(vlevels)
# }#end if
#
# #----- Loop over formats. --------------------------------------------------------#
# for (o in sequence(nout)){
# fichier = file.path(outdir,paste0(vnam,"-",suffix,".",outform[o]))
# if (outform[o] %in% "x11"){
# X11(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "quartz"){
# quartz(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "png"){
# png(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent")
# }else if (outform[o] %in% "tif"){
# tiff(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent",compression="lzw")
# }else if (outform[o] %in% "eps"){
# postscript(file=fichier,width=size$width,height=size$height
# ,pointsize=ptsz,paper=size$paper)
# }else if (outform[o] %in% "pdf"){
# pdf(file=fichier,onefile=FALSE
# ,width=size$width,height=size$height,pointsize=ptsz,paper=size$paper)
# }#end if
#
# letitre = paste0(description,"\n",lieu)
# ley = desc.unit(desc="Soil depth",unit=untab$m)
# lacle = desc.unit(desc=NULL,unit=unit)
# par(par.user)
# sombreado(x=yearaxis,y=soilaxis,z=varbuff,levels=vlevels,nlevels=vnlev
# ,color.palette=get(vcscheme)
# ,plot.title=title(main=letitre,xlab="Month",ylab=ley,cex.main=0.7)
# ,key.title=title(main=lacle,cex.main=0.8)
# ,key.log=pnlog
# ,plot.axes={axis(side=1)
# axis(side=2,at=zat,labels=znice)
# if (hovgrid){
# abline(h=zat,v=axTicks(1),col=grid.colour,lty="dotted")
# }#end if hovgrid
# }#end plot.axes
# )
#
# if (outform[o] %in% c("x11","quartz")){
# locator(n=1)
# dev.off()
# }else{
# dev.off()
# }#end if
# dummy = clean.tmp()
# } #end for outform
# }#end if plotit
# }#end for nhov
# #---------------------------------------------------------------------------------------#
#
#
#
#
#
#
# #---------------------------------------------------------------------------------------#
# # Bar plot by DBH class. (barplot_dbh folder) #
# #---------------------------------------------------------------------------------------#
# cat(" + Bar plot by DBH classes...","\n")
# pftuse = which(apply(X=szpft$nplant,MARGIN=3,FUN=sum,na.rm=TRUE) > 0.)
# pftuse = pftuse[pftuse != (npft+1)]
# npftuse = length(pftuse)
# pftname.use = pft$name [pftuse]
# pftcol.use = pft$colour[pftuse]
# my.avg = 0
# for (v in sequence(ntspftdbh)){
# #----- Load settings for this variable.----------------------------------------------#
# thisbar = tspftdbh[[v]]
# vnam = thisbar$vnam
# description = thisbar$desc
# unit = thisbar$e.unit
# stacked = thisbar$stack
# plotit = thisbar$bar.plot # && plot.ycomp
# plog = thisbar$plog
# if (plog){
# stacked = FALSE
# xylog = "y"
# }else{
# xylog = ""
# }#end if
# #------------------------------------------------------------------------------------#
#
#
# #------------------------------------------------------------------------------------#
# # Check whether to plot this
# #------------------------------------------------------------------------------------#
# if (plotit){
# cat(" - ",description,"...","\n")
#
#
# #---------------------------------------------------------------------------------#
# # Retrieve the variable, and keep only the part that is usable. #
# #---------------------------------------------------------------------------------#
# thisvnam = szpft[[vnam]]
# thisvnam = thisvnam [,,pftuse]
# thisvnam = thisvnam [,-(ndbh+1),]
#
# thisvnam[is.na(thisvnam)] = 0.
# #---------------------------------------------------------------------------------#
#
#
# #---------------------------------------------------------------------------------#
# # Find the limits for the plots. We use the same axis so it is easier to #
# # compare different times. #
# #---------------------------------------------------------------------------------#
# if (stacked){
# ylimit = c(0,max(apply(X=thisvnam,MARGIN=c(1,2),FUN=sum,na.rm=TRUE)))
# }else{
# ylimit = range(x=thisvnam,na.rm=TRUE)
# }#end if
# ylimit = pretty.xylim(u=ylimit,fracexp=0.0,is.log=plog)
# #---------------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------------#
# # Check if the directory exists. If not, create it. #
# #---------------------------------------------------------------------------------#
# barplotdir = file.path(outpref,"barplot_dbh")
# if (! file.exists(barplotdir)) dir.create(barplotdir)
# outdir = file.path(barplotdir,vnam)
# if (! file.exists(outdir)) dir.create(outdir)
# #---------------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------------#
# # Loop over all possible months. #
# #---------------------------------------------------------------------------------#
# for (y in sequence(nyears)){
#
# #----- Find which year we are plotting. ---------------------------------------#
# cyear = sprintf("%4.4i",datum$toyear[y])
# yy = as.numeric(cyear)
# #------------------------------------------------------------------------------#
#
#
# #----- Loop over output formats. ----------------------------------------------#
# for (o in sequence(nout)){
# #------ Open the plot. -----------------------------------------------------#
# fichier = file.path( outdir
# , paste0(vnam,"-",cyear,"-",suffix,".",outform[o])
# )#end file.path
# if (outform[o] %in% "x11"){
# X11(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "quartz"){
# quartz(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "png"){
# png(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth)
# }else if (outform[o] %in% "eps"){
# postscript(file=fichier,width=size$width,height=size$height
# ,pointsize=ptsz,paper=size$paper)
# }else if (outform[o] %in% "pdf"){
# pdf(file=fichier,onefile=FALSE
# ,width=size$width,height=size$height,pointsize=ptsz,paper=size$paper)
# }#end if
# #---------------------------------------------------------------------------#
#
# #----- Split plotting area into two. ---------------------------------------#
# par(par.user)
# layout(mat=rbind(2,1),heights=c(5,1))
# #---------------------------------------------------------------------------#
#
#
#
# #----- Plot the legend. ----------------------------------------------------#
# par(mar=c(0.1,4.6,0.1,2.1))
# plot.new()
# plot.window(xlim=c(0,1),ylim=c(0,1))
# legend( x = "bottom"
# , inset = 0.0
# , legend = pftname.use
# , fill = pftcol.use
# , ncol = min(3,pretty.box(n.selpft)$ncol)
# , title = expression(bold("Plant functional type"))
# , cex = cex.ptsz
# , bg = background
# , xpd = TRUE
# )
# #---------------------------------------------------------------------------#
#
#
# #------ Set up the title and axis labels. ----------------------------------#
# letitre = paste0(lieu,"\n",description," - Year : ",cyear)
# lexlab = "DBH Classes"
# leylab = desc.unit(desc=description,unit=unit)
# #---------------------------------------------------------------------------#
#
#
# #------ Correct data matrix for log plot (0s create problems, NAs not) -----#
# if (xylog == "y"){
# temp = thisvnam[y,,]
# temp[temp==0] = NA
# thisvnam[y,,] = temp
# }
# #---------------------------------------------------------------------------#
#
#
# #------- Plot all monthly means together. ------------------------------------#
# par(mar=c(4.1,4.6,4.1,2.1))
# barplot(height=t(thisvnam[y,,]),names.arg=dbhnames[1:ndbh],width=1.0,
# main=letitre,xlab=lexlab,ylab=leylab,ylim=ylimit * 1.05,
# legend.text=FALSE, beside=(! stacked),col=pftcol.use,log=xylog,
# border=grey.fg,xpd=FALSE,cex.main=cex.main,las=1)
# if (plotgrid & (! stacked)){
# xgrid=0.5+(1:ndbh)*(1+npftuse)
# abline(v=xgrid,col=grid.colour,lty="solid")
# }
# box()
# #---------------------------------------------------------------------------#
# if( v == 24 & plot.nplant.hystogram & y >= (nyears-100)){
# # dev.off()
# xylog = ""
# ylimit=c(0,140.0)
#
# mytable = thisvnam[y,,]
# mytable[is.na(mytable)] = 0.0
# my.avg=my.avg + mytable
# if (y == nyears){
#
# #multiply by ha conversion and divide by 100 (=nyears) (hence *100)
# my.avg = my.avg*100
#
# exp_lianas_dist = c(1000,117.1,53.594,28.79,15.454,
# 12.631,9.331,5.552,5.357,3.728,
# 2.816,1.188,0.754,0.573,0.5,0.5
# ,0.5)
# exp_lianas_dist_err = exp_lianas_dist * 0.3 + 2
# exp_lianas_undist = c(1000,52,31.831,20.801,10.735,
# 10.29,5.996,4.11,2.701,2.258,
# 1.814,1.369,0.5,0.5,0.5,0.5
# ,0.5)
# exp_lianas_undist_err = exp_lianas_undist * 0.4 + 2
#
# my.avg[,1] = exp_lianas_dist
# my.avg[,2] = exp_lianas_undist
#
# my.avg=my.avg[,-3]
#
# er_mat = matrix(c(exp_lianas_dist_err, exp_lianas_undist_err), nrow = 17)
# er_mat = cbind (er_mat,NA)
# er_mat [2,2] = er_mat[2,2] * 1.3
#
#
# fichier=paste(outroot,"/paracou/yearly/barplot_dbh/nplant/average.pdf",sep='')
# pdf(file=fichier,onefile=FALSE
# ,width=8.8,height=6.8,pointsize=16,paper="special")
# #----- Split plotting area into two. ---------------------------------------#
# par(par.user)
# layout(mat=rbind(2,1),heights=c(5,1))
# #---------------------------------------------------------------------------#
#
#
#
# #----- Plot the legend. ----------------------------------------------------#
# par(mar=c(0.1,4.6,0.1,2.1))
# plot.new()
# plot.window(xlim=c(0,1),ylim=c(0,1))
# legend( x = "bottom"
# , inset = 0.0
# , legend = c("Disturbed plots","Undistrurbed plots","Model")
# , fill = c("#D3D3D3","#838B8B","#1E64C8")
# , ncol = min(3,pretty.box(n.selpft)$ncol)
# , title = expression(bold("Type"))
# , cex = cex.ptsz
# , bg = background
# , xpd = TRUE
# )
# #---------------------------------------------------------------------------#
#
#
# #------ Set up the title and axis labels. ----------------------------------#
# lexlab = "DBH Classes"
# leylab = desc.unit(desc="Plant density",unit="p*l*a*n*t^phantom(1)*h*a^{-1}")
# #---------------------------------------------------------------------------#
#
# pftcol.use = c("#D3D3D3","#838B8B","#1E64C8")
#
# #------- Plot all monthly means together. ------------------------------------#
# par(mar=c(4.1,4.6,4.1,2.1))
# mynames = c("<2", "2-3", "3-4", "4-5", "5-6", "6-7", "7-8",
# "8-9", "9-10", "10-11", "11-12", " 12-13", "13-14",
# "14-15", "15-16", "16-20", ">20")
# lexlab="DBH Classes [cm]"
# barx = barplot(height=t(my.avg),names.arg=mynames,width=1.0,
# main="Diameter distribution",xlab=lexlab,ylab=leylab,ylim=ylimit,
# legend.text=FALSE, beside=(! stacked),col=pftcol.use,log="",
# border=grey.fg,xpd=FALSE,cex.main=cex.main,las=1)
# if (plotgrid & (! stacked)){
# npftuse=3
# xgrid=0.5+(1:ndbh)*(1+npftuse)
# abline(v=xgrid,col=grid.colour,lty="solid")
# }
# arrows(barx, y0=t(my.avg)-t(er_mat)/2 , y1=t(my.avg)+t(er_mat) / 2, barx, length=0)
# # arrows(barx, y1=t(my.avg)+t(er_mat), barx, height=t(my.avg), angle=90, code=1, length=0)
# box()
# #---------------------------------------------------------------------------#
#
# }
# }
#
#
#
# #---------------------------------------------------------------------------#
# # Close the device. #
# #---------------------------------------------------------------------------#
# if (outform[o] %in% c("x11","quartz")){
# locator(n=1)
# dev.off()
# }else{
# dev.off()
# }#end if
# dummy = clean.tmp()
# #---------------------------------------------------------------------------#
# } #end for outform
# #------------------------------------------------------------------------------#
# }#end for
# #---------------------------------------------------------------------------------#
# }#end if
# #------------------------------------------------------------------------------------#
# }#end for
# #---------------------------------------------------------------------------------------#
#
#---------------------------------------------------------------------------------------#
# Plot the 3-D size and age structure of various variables. (sas folder) #
#---------------------------------------------------------------------------------------#
#for (v in sequence(ntspftdbh)){
v=69
#----- Retrieve variable information from the list. ---------------------------------#
thissas = tspftdbh[[v]]
vnam = thissas$vnam
description = thissas$desc
unit = thissas$i.unit
plotit = thissas$sas
plog = thissas$plog
#----- If this variable is to be plotted, then go through this if block. ------------#
if (plotit){
cat(" + Size and age structure plot: ",description,"...","\n")
#---------------------------------------------------------------------------------#
# Check if the directory exists. If not, create it. #
#---------------------------------------------------------------------------------#
sasdir = file.path(outpref,"sas")
if (! file.exists(sasdir)) dir.create(sasdir)
outdir = file.path(sasdir,vnam)
if (! file.exists(outdir)) dir.create(outdir)
#---------------------------------------------------------------------------------#
#----- Load this list into "thislist". -------------------------------------------#
varco = cohort[[vnam]]
#---------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------#
# Loop over all times. #
#---------------------------------------------------------------------------------#
for (ww in names(cohort$age)){
#----- Find which year we are plotting. ---------------------------------------#
cmonth = substring(ww,7,8)
thisyear = substring(ww,2,5)
mm = as.numeric(cmonth)
yy = as.numeric(thisyear)
thisyear = as.character(yy - 350)
#------------------------------------------------------------------------------#
mxww = numeric()
myww = numeric()
mzww = numeric()
mpchww = numeric()
mcolww = numeric()
mcexww = numeric()
#----- Retrieve variable list, age, DBH, and PFT for this year. ---------------#
ageww = cohort$age [[ww]]
if (any(ageww <= 0,na.rm=TRUE)){
minww = min(ageww,na.rm=TRUE)
ageww = ageww - minww + 0.01
}#end if
dbhww = cohort$dbh [[ww]]
pftww = cohort$pft [[ww]]
varww = varco [[ww]]
popww = cohort$nplant[[ww]] * cohort$area[[ww]]
#------------------------------------------------------------------------------#
#------------------------------------------------------------------------------#
# We only plot the SAS figures when the polygon is not an absolute desert. #
#------------------------------------------------------------------------------#
if (any (! is.na(varww))){
#---------------------------------------------------------------------------#
# Find the range. If the user wants the range to be fixed, then use #
# the global range, otherwise, simply use the range for this year. #
#---------------------------------------------------------------------------#
if (sasfixlimits){
xlimit = pretty.xylim(u=unlist(cohort$age),fracexp=0.0,is.log=TRUE )
ylimit = pretty.xylim(u=unlist(cohort$dbh),fracexp=0.0,is.log=FALSE)
zlimit = pretty.xylim(u=unlist(varco) ,fracexp=0.0,is.log=plog )
popmin = min(unlist(cohort$nplant * cohort$area), na.rm=TRUE)
popmax = max(unlist(cohort$nplant * cohort$area), na.rm=TRUE)
}else{
xlimit = pretty.xylim(u=ageww ,fracexp=0.0,is.log=TRUE )
ylimit = pretty.xylim(u=dbhww ,fracexp=0.0,is.log=FALSE)
zlimit = pretty.xylim(u=varww ,fracexp=0.0,is.log=plog )
popmin = min(popww ,na.rm=TRUE)
popmax = max(popww ,na.rm=TRUE)
}#end if
#---------------------------------------------------------------------------#
#----- Define the scale-dependent population size. -------------------------#
cexww = cexmin + (cexmax - cexmin) * log(popww/popmin) / log(popmax/popmin)
#---------------------------------------------------------------------------#
#----- Define the floor location. ------------------------------------------#
if ((zlimit[1] > 0) != (zlimit[2] > 0)){
floor3d = 0.
}else if (zlimit[1] > 0){
floor3d = zlimit[1]
}else{
floor3d = zlimit[2]
}#end if
#---------------------------------------------------------------------------#
xlimit = c(0.1, 300)
ylimit = c(0.01, 130)
zlimit = c(0.4999, 34.999)
xlimit.stat = c(-2.5, 300)
ylimit.stat = c(0.01, 130)
zlimit.stat = c(0.4999, 34.999)
floor3d = 0.0
#----- Define the grid information for the 3-D plot. -----------------------#
xlabels = pretty.log(xlimit,n=5)
ylabels = pretty(ylimit,n=5)
zlabels = if(plog){pretty.log(zlimit,n=5)}else{pretty(zlimit,n=5)}
xat = log(xlabels)
yat = ylabels
zat = if(plog){log(zlabels)}else{zlabels}
xlimit = range(x=xat)
ylimit = range(x=yat)
zlimit = range(x=zat)
xfloor = seq(from=xlimit[1],to=xlimit[2],length.out=16)
yfloor = seq(from=ylimit[1],to=ylimit[2],length.out=16)
zfloor = matrix(floor3d,nrow=length(xfloor),ncol=length(yfloor))
#---------------------------------------------------------------------------#
#----- Expand the lines to make the lollipops. -----------------------------#
ncohnow = length(varww)
ageww = rep(ageww,each=3)
dbhww = rep(dbhww,each=3)
pftww = rep(pftww,each=3)
varww = as.vector( rbind( rep(floor3d,times=ncohnow)
, varco[[ww]]
, rep(NA,times=ncohnow)
)#end rbind
)#end as.vector
xww = log(ageww)
yww = dbhww
zww = if(plog){log(varww)}else{varww}
pchww = rep(c(NA,16,NA),times=ncohnow)
cexww = rep(cexww,each=3)
colww = pft$colour[pftww]
pftin = sort(unique(cohort$pft[[ww]]))
colleg = pft$colour[pftin]
pftleg = pft$name [pftin]
#---------------------------------------------------------------------------#
#---------------------------------------------------------------------------#
# Plot annotation. #
#---------------------------------------------------------------------------#
letitre = paste(description," - ",lieu
,"\n Time :",mlist[mm],"/",thisyear,sep=" ")
lexlab = desc.unit(desc="Gap age",unit=untab$yr)
leylab = desc.unit(desc="DBH",unit=untab$cm)
lezlab = desc.unit(desc=description,unit=unit)
#---------------------------------------------------------------------------#
#----- Loop over output formats. -------------------------------------------#
for (o in sequence(nout)){
#----- Open file. -------------------------------------------------------#
fichier = file.path( outdir
, paste0( vnam,"-",thisyear,"-",cmonth,"-",suffix
,".",outform[o]
)#end paste0
)#end file.path
if (outform[o] %in% "x11"){
X11(width=size$width,height=size$height,pointsize=ptsz)
}else if (outform[o] %in% "quartz"){
quartz(width=size$width,height=size$height,pointsize=ptsz)
}else if(outform[o] %in% "png"){
png(filename=fichier,width=size$width*depth,height=size$height*depth
,pointsize=ptsz,res=depth,bg="transparent")
}else if(outform[o] %in% "tif"){
tiff(filename=fichier,width=size$width*depth,height=size$height*depth
,pointsize=ptsz,res=depth,bg="transparent",compression="lzw")
}else if(outform[o] %in% "eps"){
postscript(file=fichier,width=size$width,height=size$height
,pointsize=ptsz,paper=size$paper)
}else if(outform[o] %in% "pdf"){
pdf(file=fichier,onefile=FALSE,width=size$width,height=size$height
,pointsize=ptsz,paper=size$paper)
}#end if
#------------------------------------------------------------------------#
#----- Split the domain into 2. -----------------------------------------#
par(par.user)
layout(mat=rbind(2,1),heights=c(5,1))
#------------------------------------------------------------------------#
#------------------------------------------------------------------------#
# Plot legend. #
#------------------------------------------------------------------------#
par(mar=c(0.1,0.1,0.1,0.1))
plot.new()
plot.window(xlim=c(0,1),ylim=c(0,1))
legend( x = "center"
, inset = 0.0
, legend = pftleg
, fill = colleg
, ncol = min(4,pretty.box(length(pftleg))$ncol)
, title = expression(bold("Plant functional type"))
, cex = cex.ptsz
, xpd = TRUE
, bty = "n"
)#end legend
#------------------------------------------------------------------------#
#------------------------------------------------------------------------#
# Plot the 3-D plot. #
#------------------------------------------------------------------------#
par(mar=c(1.1,1.1,4.1,1.1))
pout = perspx( x = xfloor
, y = yfloor
, z = zfloor
, xlim = xlimit
, ylim = ylimit
, zlim = zlimit
, theta = theta
, phi = phi
, col = gcol
, expand = expz
, ticktype = "detailed"
, border = NA
, shade = shade
, ltheta = ltheta
, main = letitre
, cex.main = 0.8*cex.ptsz
, axes = FALSE
)#end perspx
#----- Add axes. --------------------------------------------------------#
paxis3d(edge="X--",pmat=pout,at=xat,cex=0.9*cex.ptsz,labels=xlabels)
paxis3d(edge="Y--",pmat=pout,at=yat,cex=0.9*cex.ptsz,labels=ylabels)
paxis3d(edge="Z-+",pmat=pout,at=zat,cex=0.9*cex.ptsz,labels=zlabels)
mtext3d(edge="X--",pmat=pout,labels=lexlab,cex=cex.ptsz,srt=theta+90)
mtext3d(edge="Y--",pmat=pout,labels=leylab,cex=cex.ptsz,srt=theta)
mtext3d(edge="Z-+",pmat=pout,labels=lezlab,cex=cex.ptsz,srt=-75)
#------------------------------------------------------------------------#
for (i in sequence(length(xww))){
if (xww[i] %wr% xlimit.stat && yww[i] %wr% ylimit.stat){
mxww = c(mxww , xww[i])
myww = c(myww , yww[i])
mzww = c(mzww , zww[i])
mpchww = c(mpchww, pchww[i])
mcolww = c(mcolww, colww[i])
mcexww = c(mcexww, cexww[i])
}
}
if(v==69){
#----- Add the cohorts. -------------------------------------------------#
lines (trans3d(x=mxww,y=myww,z=mzww,pmat=pout),type="l",col=grey.fg,lwd=2)
points(trans3d(x=mxww,y=myww,z=mzww,pmat=pout),type="p",pch=mpchww
,col=mcolww,cex=mcexww)
#------------------------------------------------------------------------#
}
#----- Add the cohorts. -------------------------------------------------#
# lines (trans3d(x=xww,y=yww,z=zww,pmat=pout),type="l",col=grey.fg,lwd=2)
# points(trans3d(x=xww,y=yww,z=zww,pmat=pout),type="p",pch=pchww
# ,col=colww,cex=cexww)
#------------------------------------------------------------------------#
#----- Close the device. ------------------------------------------------#
if (outform[o] %in% c("x11","quartz")){
locator(n=1)
dev.off()
}else{
dev.off()
}#end if
dummy = clean.tmp()
#------------------------------------------------------------------------#
}#end for outform
#---------------------------------------------------------------------------#
}#end if is.na(varww)
#------------------------------------------------------------------------------#
}#end for nameco
#---------------------------------------------------------------------------------#
}#end if
#------------------------------------------------------------------------------------#
#}#end for npsas
#---------------------------------------------------------------------------------------#
#}#end for places
#q("no")
manfredo89/ED2io documentation built on May 21, 2019, 11:24 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#==========================================================================================#
#==========================================================================================#
# Leave these commands at the beginning. They will refresh the session. #
#------------------------------------------------------------------------------------------#
rm(list=ls())
graphics.off()
#==========================================================================================#
#==========================================================================================#
#==========================================================================================#
#==========================================================================================#
# Here is the user defined variable section. #
#------------------------------------------------------------------------------------------#
#----- Paths. -----------------------------------------------------------------------------#
args <- commandArgs(TRUE)
arg.runtype <- as.character(args[])
if (length(args) > 1){
cat("Only one directory is supported now \n")
stop()
}
if (length(args)==0) {
cat("No arguments were passed, defaulting to lianas \n")
arg.runtype="lianas"
}
run.type = arg.runtype
here = getwd() # Current directory.
there = paste(here,"/../",run.type[1],sep='') # Directory where analyses/history are
srcdir = here # Source directory.
outroot = paste(there,"/figures",sep='') # Directory for figures
cat("here:",here,"\n")
cat("there:",there,"\n")
cat("srcdir:",srcdir,"\n")
cat("outroot:",outroot,"\n")
#------------------------------------------------------------------------------------------#
#----- Time options. ----------------------------------------------------------------------#
monthbeg = 01 # First month to use
yearbeg = 2000 # First year to consider
yearend = 2050 # Maximum year to consider
reload.data = T # Should I reload partially loaded data?
#sasmonth.short = c(2,5,8,11) # Months for SAS plots (short runs)
sasmonth.short = c(2,4,6,8,10,12) # Months for SAS plots (short runs)
sasmonth.long = c(2,4,6,8,10,12) # Months for SAS plots (short runs)
#sasmonth.long = 5 # Months for SAS plots (long runs)
#nyears.long = 800 # Runs longer than this are considered long runs.
#------------------------------------------------------------------------------------------#
#----- Name of the simulations. -----------------------------------------------------------#
myplaces = c("paracou")
#------------------------------------------------------------------------------------------#
#----- Plot options. ----------------------------------------------------------------------#
outform = c("pdf") # Formats for output file. Supported formats are:
# - "X11" - for printing on screen
# - "quartz" - for printing on Mac OS screen
# - "eps" - for postscript printing
# - "png" - for PNG printing
# - "tif" - for TIFF printing
# - "pdf" - for PDF printing
depth = 96 # PNG resolution, in pixels per inch
paper = "letter" # Paper size, to define the plot shape
ptsz = 16 # Font size.
lwidth = 2.5 # Line width
plotgrid = TRUE # Should I plot the grid in the background?
sasfixlimits = FALSE # Use a fixed scale for size and age-structure
# plots? (FALSE will set a suitable scale for
# each plot)
fcgrid = TRUE # Include a grid on the filled contour plots?
ncolshov = 200 # Target number of colours for Hovmoller diagrams.
hovgrid = TRUE # Include a grid on the Hovmoller plots?
legwhere = "topleft" # Where should I place the legend?
inset = 0.01 # Inset between legend and edge of plot region.
scalleg = 0.40 # Expand y limits by this relative amount to fit
# the legend
cex.main = 0.8 # Scale coefficient for the title
theta = 315. # Azimuth for perspective projection
phi = 30. # Vertical angle for perspective projection
ltheta = -210. # Azimuth angle for light
shade = 0.125 # Shade intensity
expz = 0.5 # Expansion factor for Z axis
cexmin = 0.5 # Minimum "head" size of the lollipop
cexmax = 3.0 # Maximum "head" size of the lollipop
ylnudge = 0.05 # Nudging factor for ylimit
ptype = "l" # Type of plot
ptyped = "p" # Type of plot
ptypeb = "o" # Type of plot
drought.mark = FALSE # Put a background to highlight droughts?
drought.yeara = 1605 # First year that has drought
drought.yearz = 1609 # Last year that has drought
months.drought = c(12,1,2,3) # Months with drought
ibackground = 0 # Background settings (check load_everything.r)
plot.nplant.hystogram = T # Plot custom dbh distribution (works but bad implementation, correct or remove)
options(bitmapType='cairo')
#------------------------------------------------------------------------------------------#
#------ Miscellaneous settings. -----------------------------------------------------------#
slz.min = -5.0 # The deepest depth that trees access water.
idbh.type = 3 # Type of DBH class
# 1 -- Every 10 cm until 100cm; > 100cm
# 2 -- 0-10; 10-20; 20-35; 35-50; 50-70; > 70 (cm)
# 3 -- 0-10; 10-35; 35-55; > 55 (cm)
klight = 0.8 # Weighting factor for maximum carbon balance
corr.growth.storage = 1.0 # Correction factor to be applied to growth and
# storage respiration
iallom = 3 # Allometry to use
#------------------------------------------------------------------------------------------#
#==========================================================================================#
#==========================================================================================#
#==========================================================================================#
#==========================================================================================#
#==========================================================================================#
#==========================================================================================#
# NO NEED TO CHANGE ANYTHING BEYOND THIS POINT UNLESS YOU ARE DEVELOPING THE CODE... #
#==========================================================================================#
#==========================================================================================#
#==========================================================================================#
#==========================================================================================#
#==========================================================================================#
#==========================================================================================#
analy.path = paste(there,"analy",sep='/') # analysis folder
com.list = list.files(analy.path, pattern = "*-Q-*")
if(length(run.type) > 1) com.list = com.list[duplicated(com.list)]
if(length(com.list) < 24){ stop("ERROR: two few files in analysis folder") }
#------------------------------------------------------------------------------------------#
# Keep only full years #
#------------------------------------------------------------------------------------------#
# first year of available data in the folder
if (is.na(yearbeg)) yearbeg = min(sub('\\-.*','',sub("^.*?Q-","",com.list)))
# last year if available data in the folder
if (is.na(yearend)) yearend = max(sub('\\-.*','',sub("^.*?Q-","",com.list)))
yearbeg = as.integer(yearbeg) # convert yeara to factor
yearend = as.integer(yearend) # convert yearz to factor
#----- Load some packages and scripts. ----------------------------------------------------#
source(file.path(srcdir,"load.everything.r"))
#------------------------------------------------------------------------------------------#
#----- Set how many formats we must output. -----------------------------------------------#
outform = tolower(outform)
nout = length (outform)
#------------------------------------------------------------------------------------------#
#----- Avoid unecessary and extremely annoying beeps. -------------------------------------#
options(locatorBell=FALSE)
#------------------------------------------------------------------------------------------#
#----- Load observations. -----------------------------------------------------------------#
#obsrfile = file.path(srcdir,"LBA_MIP.v8.RData")
#load(file=obsrfile)
#----- Define plot window size ------------------------------------------------------------#
size = plotsize(proje=FALSE,paper=paper)
#------------------------------------------------------------------------------------------#
#---- Create the main output directory in case there is none. -----------------------------#
if (! file.exists(outroot)) dir.create(outroot)
#------------------------------------------------------------------------------------------#
#------------------------------------------------------------------------------------------#
# Big place loop starts here... #
#------------------------------------------------------------------------------------------#
#for (place in myplaces){
place="paracou"
#----- Retrieve default information about this place and set up some variables. --------#
thispoi = locations(where=place,here=there,yearbeg=yearbeg,yearend=yearend
,monthbeg=monthbeg)
inpref = thispoi$pathin
outmain = file.path(outroot,place)
outpref = file.path(outmain,"yearly")
lieu = thispoi$lieu
iata = thispoi$iata
suffix = thispoi$iata
yeara = thispoi$yeara
yearz = thispoi$yearz
meszz = thispoi$monz
#---------------------------------------------------------------------------------------#
# Make sure we only deal with full years. #
#---------------------------------------------------------------------------------------#
if (monthbeg > 1) yeara = yeara + 1
if (meszz < 12) yearz = yearz - 1
monthbeg = 1
meszz = 12
if (yeara > yearz){
cat(" - Yeara: ",yeara,"\n")
cat(" - Yearz: ",yearz,"\n")
cat(" - Prefix: ",inpref,"\n")
cat(" - Invalid years, will not process data...","\n")
q("no")
}#end if
#---------------------------------------------------------------------------------------#
#----- Create the directories in case they don't exist. --------------------------------#
if (! file.exists(outmain)) dir.create(outmain)
if (! file.exists(outpref)) dir.create(outpref)
#---------------------------------------------------------------------------------------#
#----- Decide how frequently the cohort-level variables should be saved. ---------------#
if ((yearend - yearbeg + 1) <= nyears.long){
sasmonth = sasmonth.short
plot.ycomp = TRUE
}else{
sasmonth = sasmonth.long
# plot.ycomp = FALSE
plot.ycomp = TRUE
}#end if
#---------------------------------------------------------------------------------------#
#----- Print a banner to entretain the user. -------------------------------------------#
cat(" + Post-processing output from ",lieu,"...","\n")
#---------------------------------------------------------------------------------------#
# Flush all variables that will hold the data. #
#---------------------------------------------------------------------------------------#
ntimes = (yearz-yeara-1)*12+meszz+(12-monthbeg+1)
nyears = yearz-yeara+1
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#
# Make the RData file name, then we check whether we must read the files again #
# or use the stored RData. Notice that the path is the same for plot_ycomp.r and #
# plot_monthly, so you don't need to read in the data twice. #
#---------------------------------------------------------------------------------------#
path.data = file.path(there,"rdata_month")
if (! file.exists(path.data)) dir.create(path.data)
ed22.rdata = file.path(path.data,paste("M",place,"RData",sep="."))
ed22.status = file.path(path.data,paste("Mstatus_",place,".txt",sep=""))
if (reload.data && file.exists(ed22.rdata)){
#----- Load the modelled dataset. ---------------------------------------------------#
cat(" - Loading previous session...","\n")
load(ed22.rdata)
tresume = datum$ntimes + 1
if (ntimes > datum$ntimes){
datum = update.monthly( new.ntimes = ntimes
, old.datum = datum
, montha = monthbeg
, yeara = yeara
, inpref = inpref
, slz.min = slz.min
)#end update.monthly
}#end if
#------------------------------------------------------------------------------------#
}else{
cat(" - Starting new session...","\n")
tresume = 1
datum = create.monthly( ntimes = ntimes
, montha = monthbeg
, yeara = yeara
, inpref = inpref
, slz.min = slz.min
)#end create.monthly
}#end if
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#
# Check whether we have anything to update. #
#---------------------------------------------------------------------------------------#
complete = tresume > ntimes
#---------------------------------------------------------------------------------------#
#----- Copy some dimensions to scalars. ------------------------------------------------#
nzg = datum$nzg
nzs = datum$nzs
ndcycle = datum$ndcycle
isoilflg = datum$isoilflg
slz = datum$slz
slxsand = datum$slxsand
slxclay = datum$slxclay
ntext = datum$ntext
soil.prop = datum$soil.prop
dslz = datum$dslz
soil.depth = datum$soil.depth
soil.dry = datum$soil.dry
soil.poro = datum$soil.poro
ka = datum$ka
kz = datum$kz
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#
# Loop over all times in case there is anything new to be read. #
#---------------------------------------------------------------------------------------#
if (! complete){
#------------------------------------------------------------------------------------#
# This function will read the files. #
#------------------------------------------------------------------------------------#
datum = read.q.files(datum=datum,ntimes=ntimes,tresume=tresume,sasmonth=sasmonth)
#------------------------------------------------------------------------------------#
#------ Save the data to the R object. ----------------------------------------------#
cat(" + Saving data to ",basename(ed22.rdata),"...","\n")
save(datum,file=ed22.rdata)
#------------------------------------------------------------------------------------#
}#end if (! complete)
#---------------------------------------------------------------------------------------#
#----- Update status file with latest data converted into R. ---------------------------#
#latest = paste(datum$year[ntimes],datum$month[ntimes],sep=" ")
#dummy = write(x=latest,file=ed22.status,append=FALSE)
#---------------------------------------------------------------------------------------#
#----- Make some shorter versions of some variables. -----------------------------------#
mfac = datum$month
yfac = datum$year
emean = datum$emean
emsqu = datum$emsqu
qmean = datum$qmean
qmsqu = datum$qmsqu
szpft = datum$szpft
lu = datum$lu
patch = datum$patch
cohort = datum$cohort
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#
# Consolidate the yearly means for the long-term dynamics (the PFT and DBH/PFT #
# stuff). #
#---------------------------------------------------------------------------------------#
cat (" - Finding the annual statistics for multi-dimensional variables...","\n")
cat (" * Aggregating the annual mean of PFT-DBH variables...","\n")
for (vname in names(szpft)){
szpft[[vname]] = qapply(X=szpft[[vname]],INDEX=yfac,DIM=1,FUN=mean,na.rm=TRUE)
}#end for
#----- LU arrays. The "+1" column contains the total. --------------------------------#
cat (" * Aggregating the annual mean of LU variables...","\n")
for (vname in names(lu)){
lu [[vname]] = qapply(X=lu [[vname]],INDEX=yfac,DIM=1,FUN=mean,na.rm=TRUE)
}#end for
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#
# Here we find the monthly means for month, then compute the standard deviation. #
#---------------------------------------------------------------------------------------#
cat (" - Finding the monthly and annual means...","\n")
cat (" * Aggregating the monthly mean and standard deviation...","\n")
mmean = list()
msdev = list()
ymean = list()
ysdev = list()
for (vname in names(emean)){
if (vname %in% c("soil.temp","soil.water","soil.mstpot","soil.extracted")){
mmean[[vname]] = qapply(X=emean[[vname]], INDEX=mfac, DIM=1, FUN=mean, na.rm=TRUE)
msdev[[vname]] = qapply(X=emean[[vname]], INDEX=mfac, DIM=1, FUN=sd , na.rm=TRUE)
ymean[[vname]] = qapply(X=emean[[vname]], INDEX=yfac, DIM=1, FUN=mean, na.rm=TRUE)
ysdev[[vname]] = qapply(X=emean[[vname]], INDEX=yfac, DIM=1, FUN=sd , na.rm=TRUE)
}else if (vname %in% c("rain","runoff","intercepted","wshed")){
mmean[[vname]] = tapply(X=emean[[vname]], INDEX=mfac, FUN=mean, na.rm=TRUE)
msdev[[vname]] = tapply(X=emean[[vname]], INDEX=mfac, FUN=sd , na.rm=TRUE)
ymean[[vname]] = tapply(X=emean[[vname]], INDEX=yfac, FUN=sum , na.rm=TRUE)
ysdev[[vname]] = tapply(X=emean[[vname]], INDEX=yfac, FUN=sd , na.rm=TRUE)
}else{
mmean[[vname]] = tapply(X=emean[[vname]], INDEX=mfac, FUN=mean, na.rm=TRUE)
msdev[[vname]] = tapply(X=emean[[vname]], INDEX=mfac, FUN=sd , na.rm=TRUE)
ymean[[vname]] = tapply(X=emean[[vname]], INDEX=yfac, FUN=mean, na.rm=TRUE)
ysdev[[vname]] = tapply(X=emean[[vname]], INDEX=yfac, FUN=sd , na.rm=TRUE)
}#end if
#------------------------------------------------------------------------------------#
#----- Fix the bad data. ------------------------------------------------------------#
bad.mmean = ! is.finite(mmean[[vname]])
bad.msdev = ! is.finite(msdev[[vname]])
bad.ymean = ! is.finite(ymean[[vname]])
bad.ysdev = ! is.finite(ysdev[[vname]])
mmean[[vname]][bad.mmean] = NA
msdev[[vname]][bad.msdev] = 0.
ymean[[vname]][bad.ymean] = NA
ysdev[[vname]][bad.ysdev] = 0.
#------------------------------------------------------------------------------------#
}#end for
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#
# Here we find the Mean diurnal cycle for each month, then compute the standard #
# deviation. #
#---------------------------------------------------------------------------------------#
cat (" - Aggregating the annual mean and std. dev. of the diurnal cycle...","\n")
umean = list()
usdev = list()
for (vname in names(qmean)){
umean[[vname]] = qapply(qmean[[vname]],INDEX=yfac,DIM=1,FUN=mean,na.rm=TRUE)
usdev[[vname]] = qapply(qmean[[vname]],INDEX=yfac,DIM=1,FUN=sd ,na.rm=TRUE)
bad.umean = ! is.finite(umean[[vname]])
bad.usdev = ! is.finite(usdev[[vname]])
umean[[vname]][bad.umean] = NA
usdev[[vname]][bad.usdev] = 0.
}#end for
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#
# Remove all elements of the DBH/PFT class that do not have a single valid cohort #
# at any given time. #
#---------------------------------------------------------------------------------------#
empty = is.na(szpft$nplant) | szpft$nplant == 0
for (vname in names(szpft)) szpft[[vname]][empty] = NA
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#
# Convert mortality and recruitment so it is scaled between 0 and 100%. #
#---------------------------------------------------------------------------------------#
szpft$mort = 100. * (1.0 - exp(- szpft$mort ) )
szpft$dimort = 100. * (1.0 - exp(- szpft$dimort ) )
szpft$ncbmort = 100. * (1.0 - exp(- szpft$ncbmort ) )
szpft$recrpft = 100. * ( exp( szpft$recr ) - 1.0)
szpft$agb.mort = 100. * (1.0 - exp(- szpft$agb.mort ) )
szpft$agb.dimort = 100. * (1.0 - exp(- szpft$agb.dimort ) )
szpft$agb.ncbmort = 100. * (1.0 - exp(- szpft$agb.ncbmort ) )
szpft$agb.recrpft = 100. * ( exp( szpft$agb.recr ) - 1.0)
szpft$bsa.mort = 100. * (1.0 - exp(- szpft$bsa.mort ) )
szpft$bsa.dimort = 100. * (1.0 - exp(- szpft$bsa.dimort ) )
szpft$bsa.ncbmort = 100. * (1.0 - exp(- szpft$bsa.ncbmort ) )
szpft$bsa.recrpft = 100. * ( exp( szpft$bsa.recr ) - 1.0)
#---------------------------------------------------------------------------------------#
#----- Find which PFTs, land uses and transitions we need to consider ------------------#
pftave = apply( X = szpft$agb[,ndbh+1,]
, MARGIN = 2
, FUN = mean
, na.rm = TRUE
)#end apply
luave = apply( X = lu$agb
, MARGIN = 2
, FUN = mean
, na.rm = TRUE
)#end apply
distave = apply(X=lu$dist,MARGIN=c(2,3),FUN=mean)
selpft = is.finite(pftave ) & pftave > 0.
sellu = is.finite(luave ) & luave > 0.
seldist = is.finite(distave) & distave > 0.
n.selpft = sum(selpft )
n.sellu = sum(sellu )
n.seldist = sum(seldist)
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#
# Define a suitable scale for diurnal cycle... #
#---------------------------------------------------------------------------------------#
thisday = seq(from=0,to=ndcycle,by=1) * 24 / ndcycle
uplot = list()
uplot$levels = c(0,4,8,12,16,20,24)
uplot$n = 7
uplot$scale = "hours"
uplot$padj = rep(0,times=uplot$n)
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#
# Define a suitable scale for soil profile layers... #
#---------------------------------------------------------------------------------------#
znice = -pretty.log(-slz,n=8)
znice = sort(c(znice,slz[1],slz[nzg]))
sel = znice >= slz[1] & znice <= slz[nzg]
znice = znice[sel]
zat = -log(-znice)
nznice = length(znice)
znice = sprintf("%.2f",znice)
#---------------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------------#
# Define a suitable scale for monthly means... #
#---------------------------------------------------------------------------------------#
montmont = seq(from=1,to=12,by=1)
mplot = list()
mplot$levels = montmont
mplot$labels = capwords(mon2mmm(montmont))
mplot$n = 12
mplot$scale = "months"
mplot$padj = rep(0,times=mplot$n)
#---------------------------------------------------------------------------------------#
#=======================================================================================#
#=======================================================================================#
#=======================================================================================#
# Plotting section begins here... #
#---------------------------------------------------------------------------------------#
cat (" - Plotting figures...","\n")
#---------------------------------------------------------------------------------------#
#
# #---------------------------------------------------------------------------------------#
# # Time series by PFT. (tspft folder) #
# #---------------------------------------------------------------------------------------#
# for (v in sequence(ntspftdbh)){
# thistspft = tspftdbh[[v]]
# vnam = thistspft$vnam
# description = thistspft$desc
# unit = thistspft$e.unit
# plog = thistspft$plog
# plotit = thistspft$pft
#
# #----- Check whether the user wants to have this variable plotted. ------------------#
# if (plotit && any(selpft)){
#
# #---------------------------------------------------------------------------------#
# # Check whether the time series directory exists. If not, create it. #
# #---------------------------------------------------------------------------------#
# outdir = file.path(outpref,"tspft")
# if (! file.exists(outdir)) dir.create(outdir)
# cat(" +",description,"time series for all PFTs...","\n")
#
# #----- Load variable -------------------------------------------------------------#
# if (vnam %in% names(szpft)){
# thisvar = szpft[[vnam]][,ndbh+1,]
# if (plog){
# #----- Eliminate non-positive values in case it is a log plot. -------------#
# badlog = is.finite(thisvar) & thisvar <= 0
# thisvar[badlog] = NA
# }#end if
# }else{
# thisvar = matrix(NA,ncol=npft+1,nrow=nyears)
# }#end if
# #---------------------------------------------------------------------------------#
#
#
#
# #----- Loop over output formats. -------------------------------------------------#
# for (o in sequence(nout)){
# fichier = file.path(outdir,paste0(vnam,"-",suffix,".",outform[o]))
# if(outform[o] %in% "x11"){
# X11(width=size$width,height=size$height,pointsize=ptsz)
# }else if(outform[o] %in% "quartz"){
# quartz(width=size$width,height=size$height,pointsize=ptsz)
# }else if(outform[o] %in% "png"){
# png(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent")
# }else if(outform[o] %in% "tif"){
# tiff(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent",compression="lzw")
# }else if(outform[o] %in% "eps"){
# postscript(file=fichier,width=size$width,height=size$height
# ,pointsize=ptsz,paper=size$paper)
# }else if(outform[o] %in% "pdf"){
# pdf(file=fichier,onefile=FALSE
# ,width=size$width,height=size$height,pointsize=ptsz,paper=size$paper)
# }#end if
#
#
# #------------------------------------------------------------------------------#
# # Find the limit, make some room for the legend, and in case the field is #
# # a constant, nudge the limits so the plot command will not complain. #
# #------------------------------------------------------------------------------#
# xlimit = pretty.xylim(u = datum$toyear ,fracexp=0.0,is.log=FALSE)
# #manfredo 574 instead of ylimit = pretty.xylim(u = thisvar[,selpft],fracexp=0.0,is.log=plog )
# ylimit = pretty.xylim(u = thisvar[selpft],fracexp=0.2,is.log=plog )
# if (plog){
# xylog = "y"
# ydrought = c( exp(sqrt(ylimit[1]^3/ylimit[2]))
# , exp(sqrt(ylimit[2]^3/ylimit[1]))
# )#end c
# }else{
# xylog = ""
# ydrought = c( ylimit[1] - 0.5 * diff(ylimit),ylimit[2] + 0.5 * diff(ylimit) )
# }#end if
# #------------------------------------------------------------------------------#
#
#
# #----- Plot settings. ---------------------------------------------------------#
# letitre = paste(description,lieu,sep=" - ")
# ley = desc.unit(desc=description,unit=unit)
# cols = pft$colour[selpft]
# legs = pft$name [selpft]
# #------------------------------------------------------------------------------#
#
#
# #------------------------------------------------------------------------------#
# # Split the plot into two windows. #
# #------------------------------------------------------------------------------#
# par(par.user)
# layout(mat=rbind(2,1),heights=c(5,1))
# #------------------------------------------------------------------------------#
#
#
#
# #------------------------------------------------------------------------------#
# # First plot: legend. #
# #------------------------------------------------------------------------------#
# par(mar=c(0.1,4.6,0.1,2.1))
# plot.new()
# plot.window(xlim=c(0,1),ylim=c(0,1))
# legend( x = "bottom"
# , inset = 0.0
# , legend = legs
# , col = cols
# , lwd = lwidth
# , ncol = min(pretty.box(n.selpft)$ncol,3)
# , title = expression(bold("Plant Functional Type"))
# , xpd = TRUE
# , bty = "n"
# )#end legend
# #------------------------------------------------------------------------------#
#
#
#
# #------------------------------------------------------------------------------#
# # Main plot. #
# #------------------------------------------------------------------------------#
# par(mar=c(4.1,4.6,4.1,2.1))
# plot.new()
# plot.window(xlim=xlimit,ylim=ylimit,log=xylog)
# axis(side=1)
# axis(side=2,las=1)
# box()
# #title(main=letitre,xlab="Year",ylab=ley,cex.main=0.7,log=xylog)
# title(main=letitre,xlab="Year",ylab=ley,cex.main=0.7)
# if (drought.mark){
# for (n in sequence(ndrought)){
# rect(xleft = drought[[n]][1],ybottom = ydrought[1]
# ,xright = drought[[n]][2],ytop = ydrought[2]
# ,col = grid.colour,border=NA)
# }#end for
# }#end if
# #----- Plot grid. -------------------------------------------------------------#
# if (plotgrid){
# abline(v=axTicks(side=1),h=axTicks(side=2),col=grid.colour,lty="solid")
# }#end if
# #----- Plot lines. ------------------------------------------------------------#
# for (n in sequence(npft+1)){
# if (selpft[n]){
# lines(datum$toyear,thisvar[,n],type="l",col=pft$colour[n],lwd=lwidth)
# }#end if
# }#end for
# #------------------------------------------------------------------------------#
#
#
# #----- Close the device. ------------------------------------------------------#
# if (outform[o] %in% c("x11","quartz")){
# locator(n=1)
# dev.off()
# }else{
# dev.off()
# }#end if
# dummy=clean.tmp()
# #------------------------------------------------------------------------------#
# } #end for outform
# }#end if (tseragbpft)
# } #end for tseries
# #---------------------------------------------------------------------------------------#
#
#
#
#
# #---------------------------------------------------------------------------------------#
# # Time series by DBH, by PFT. (tsdbh folder) #
# #---------------------------------------------------------------------------------------#
# #----- Find the PFTs to plot. ----------------------------------------------------------#
# pftuse = which(apply(X=is.na(szpft$nplant),MARGIN=3,FUN=sum,na.rm=TRUE) == 0.)
# pftuse = pftuse[pftuse != (npft+1)]
# for (v in sequence(ntspftdbh)){
# thistspftdbh = tspftdbh[[v]]
# vnam = thistspftdbh$vnam
# description = thistspftdbh$desc
# unit = thistspftdbh$e.unit
# plog = thistspftdbh$plog
# plotit = thistspftdbh$pftdbh
#
# #----- Load variable ----------------------------------------------------------------#
# if (vnam %in% names(szpft)){
# thisvar = szpft[[vnam]]
# if (plog){
# xylog="y"
# badlog = is.finite(thisvar) & thisvar <= 0
# thisvar[badlog] = NA
# }else{
# xylog=""
# }#end if
# }else{
# thisvar = array(NA,dim=c(nyears,ndbh+1,npft+1))
# }#end if
# #----- Check whether the user wants to have this variable plotted. ------------------#
# if (plotit && length(pftuse) > 0 && any(is.finite(thisvar))){
#
# #---------------------------------------------------------------------------------#
# # Check whether the time series directory exists. If not, create it. #
# #---------------------------------------------------------------------------------#
# outdir = file.path(outpref,"tsdbh")
# if (! file.exists(outdir)) dir.create(outdir)
# outvar = file.path(outdir,vnam)
# if (! file.exists(outvar)) dir.create(outvar)
# #---------------------------------------------------------------------------------#
#
# cat(" +",description,"time series for DBH class...","\n")
#
#
# #---------------------------------------------------------------------------------#
# # Find the limit, make some room for the legend, and in case the field is a #
# # constant, nudge the limits so the plot command will not complain. #
# #---------------------------------------------------------------------------------#
# xlimit = pretty.xylim(u=datum$toyear ,fracexp=0.0,is.log=FALSE)
# ylimit = pretty.xylim(u=thisvar[,,pftuse],fracexp=0.0,is.log=plog)
# if (plog){
# xylog = "y"
# ydrought = c( exp(sqrt(ylimit[1]^3/ylimit[2]))
# , exp(sqrt(ylimit[2]^3/ylimit[1]))
# )#end c
# }else{
# xylog = ""
# ydrought = c( ylimit[1] - 0.5 * diff(ylimit),ylimit[2] + 0.5 * diff(ylimit) )
# }#end if
# #---------------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------------#
# # Loop over plant functional types. #
# #---------------------------------------------------------------------------------#
# for (p in pftuse){
# pftlab = paste0("pft-",sprintf("%2.2i",p))
# cat(" - ",pft$name[p],"\n")
#
#
# #----- Loop over output formats. ----------------------------------------------#
# for (o in sequence(nout)){
# #----- Open file. ----------------------------------------------------------#
# fichier = file.path( outvar
# , paste0(vnam,"-",pftlab,"-",suffix,".",outform[o])
# )#end file.path
# if (outform[o] %in% "x11"){
# X11(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "quartz"){
# quartz(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "png"){
# png(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent")
# }else if (outform[o] %in% "tif"){
# tiff(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent",compression="lzw")
# }else if (outform[o] %in% "eps"){
# postscript(file=fichier,width=size$width,height=size$height
# ,pointsize=ptsz,paper=size$paper)
# }else if (outform[o] %in% "pdf"){
# pdf(file=fichier,onefile=FALSE
# ,width=size$width,height=size$height,pointsize=ptsz,paper=size$paper)
# }#end if
# #---------------------------------------------------------------------------#
#
#
#
# #----- Plot annotation. ---------------------------------------------------#
# letitre = paste(description,pft$name[p],lieu,sep=" - ")
# ley = desc.unit(desc=description,unit=unit)
# #---------------------------------------------------------------------------#
#
#
# #---------------------------------------------------------------------------#
# # Split the plot into two windows. #
# #---------------------------------------------------------------------------#
# par(par.user)
# layout(mat=rbind(2,1),heights=c(5,1))
# #---------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------#
# # First plot: legend. #
# #---------------------------------------------------------------------------#
# par(mar=c(0.1,4.6,0.1,2.1))
# plot.new()
# plot.window(xlim=c(0,1),ylim=c(0,1))
# legend( x = "bottom"
# , inset = 0.0
# , bg = background
# , legend = dbhnames
# , col = dbhcols
# , ncol = min(pretty.box(ndbh+1)$ncol,3)
# , title = expression(bold("DBH class"))
# , lwd = lwidth
# , xpd = TRUE
# , bty = "n"
# )#end legend
# #---------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------#
# # Main plot. #
# #---------------------------------------------------------------------------#
# par(mar=c(4.1,4.6,4.1,2.1))
# plot.new()
# plot.window(xlim=xlimit,ylim=ylimit,log=xylog)
# axis(side=1)
# axis(side=2,las=1)
# box()
# title(main=letitre,xlab="Year",ylab=ley,cex.main=0.7)
# if (drought.mark){
# for (n in sequence(ndrought)){
# rect(xleft = drought[[n]][1],ybottom = ydrought[1]
# ,xright = drought[[n]][2],ytop = ydrought[2]
# ,col = grid.colour,border=NA)
# }#end for
# }#end if
# #----- Plot grid. ----------------------------------------------------------#
# if (plotgrid){
# abline(v=axTicks(side=1),h=axTicks(side=2),col=grid.colour,lty="solid")
# }#end if
# #----- Plot lines. ---------------------------------------------------------#
# for (d in seq(from=1,to=ndbh+1,by=1)){
# lines(datum$toyear,thisvar[,d,p],type="l",col=dbhcols[d],lwd=lwidth)
# }#end for
# #---------------------------------------------------------------------------#
#
#
#
# #----- Close the device. ---------------------------------------------------#
# if (outform[o] %in% c("x11","quartz")){
# locator(n=1)
# dev.off()
# }else{
# dev.off()
# }#end if
# dummy=clean.tmp()
# #---------------------------------------------------------------------------#
# }#end for outform
# #------------------------------------------------------------------------------#
# }#end for (p in pftuse)
# #---------------------------------------------------------------------------------#
# }#end if (tseragbpft)
# #------------------------------------------------------------------------------------#
# } #end for tseries
# #---------------------------------------------------------------------------------------#
#
#
#
#
# #---------------------------------------------------------------------------------------#
# # Plot the comparison between observations and model. (ycomp folder) #
# #---------------------------------------------------------------------------------------#
# cat(" + Year-by-year comparisons of monthly means...","\n")
# for (cc in sequence(ncompmodel)){
#
# #----- Retrieve variable information from the list. ---------------------------------#
# compnow = compmodel[[cc]]
# vname = compnow$vnam
# description = compnow$desc
# unit = compnow$unit
# plotsd = compnow$plotsd
# lcolours = compnow$colour
# errcolours = compnow$errcol
# angle = compnow$angle
# dens = compnow$dens
# llwd = compnow$lwd
# shwd = compnow$shwd
# ltype = compnow$type
# plog = compnow$plog
# legpos = compnow$legpos
# plotit = compnow$mmean
#
# plotit = ( plotit && vname %in% names(emean) && vname %in% names(mmean)
# && plot.ycomp )
#
# if (plotit){
# #---------------------------------------------------------------------------------#
# # Copy the observations to a scratch variable. #
# #---------------------------------------------------------------------------------#
# thisvar = emean [[vname]]
# thismean = mmean [[vname]]
# if (length(msdev[[vname]]) == 0){
# thissdev = 0. * thismean
# }else{
# thissdev = msdev[[vname]]
# }#end if
# mod.x = montmont
# mod.ylow = thismean - thissdev
# mod.yhigh = thismean + thissdev
# mod.x.poly = c(mod.x,rev(mod.x))
# mod.y.poly = c(mod.ylow,rev(mod.yhigh))
# mod.keep = is.finite(mod.y.poly)
# mod.x.poly = mod.x.poly[mod.keep]
# mod.y.poly = mod.y.poly[mod.keep]
# #---------------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------------#
# # Check whether the time series directory exists. If not, create it. #
# #---------------------------------------------------------------------------------#
# outdir = file.path(outpref,"ycomp")
# outvar = file.path(outdir,vname)
# if (! file.exists(outdir)) dir.create(outdir)
# if (! file.exists(outvar)) dir.create(outvar)
# cat(" - ",description,"comparison...","\n")
# #---------------------------------------------------------------------------------#
#
#
#
# #----- Find the plot range. ------------------------------------------------------#
# if (plotsd){
# ylimit = range(c(mod.ylow,mod.yhigh,thisvar),na.rm=TRUE)
# }else{
# ylimit = range(thisvar,na.rm=TRUE)
# }#end if
# #----- Expand the upper range in so the legend doesn't hide things. --------------#
# ylimit = pretty.xylim(u=ylimit,fracexp=0.0,is.log=FALSE)
# #---------------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------------#
# # Loop over all years, and make one plot per year. #
# #---------------------------------------------------------------------------------#
# for (y in sequence(nyears)){
# #----- Retrieve the year and the variable for this year. ----------------------#
# year.now = datum$toyear[y]
# cyear = sprintf("%4.4i",year.now)
# var.year = thisvar[yfac == year.now]
# #------------------------------------------------------------------------------#
#
#
#
# #----- Load variable ----------------------------------------------------------#
# letitre = paste0(description," - ",lieu,"\n","Monthly mean - ",cyear)
# ley = desc.unit(desc=description,unit=unit)
# #------------------------------------------------------------------------------#
#
#
# #----- Loop over formats. -----------------------------------------------------#
# for (o in sequence(nout)){
# fichier = file.path(outvar,paste0(vname,"-",cyear,".",outform[o]))
# if (outform[o] %in% "x11"){
# X11(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "quartz"){
# quartz(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "png"){
# png(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent")
# }else if (outform[o] %in% "tif"){
# tiff(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent",compression="lzw")
# }else if (outform[o] %in% "eps"){
# postscript(file=fichier,width=size$width,height=size$height
# ,pointsize=ptsz,paper=size$paper)
# }else if (outform[o] %in% "pdf"){
# pdf(file=fichier,onefile=FALSE
# ,width=size$width,height=size$height,pointsize=ptsz,paper=size$paper)
# }#end if
#
#
#
# #----- Split plot into two windows. ----------------------------------------#
# par(par.user)
# layout(mat=rbind(2,1),heights=c(5,1))
# #---------------------------------------------------------------------------#
#
#
#
# #------ First plot: the legend. --------------------------------------------#
# par(mar=c(0.1,4.6,0.1,2.1))
# plot.new()
# plot.window(xlim=c(0,1),ylim=c(0,1))
# if (plotsd){
# legend( x = "bottom"
# , inset = 0.0
# , legend = c(cyear,paste0("Mean: ",yeara,"-",yearz))
# , fill = errcolours
# , angle = angle
# , density = dens
# , lwd = llwd
# , col = lcolours
# , bg = background
# , title = expression(bold("Shaded areas = 1 SD"))
# , cex = cex.ptsz
# , pch = 16
# , xpd = TRUE
# , bty = "n"
# )#end legend
# }else{
# legend( x = "bottom"
# , inset = 0.0
# , legend = c(cyear,paste0("Mean: ",yeara,"-",yearz))
# , col = lcolours
# , lwd = llwd
# , cex = cex.ptsz
# , pch = 16
# , xpd = TRUE
# , bty = "n"
# )#end legend
# }#end if
# #---------------------------------------------------------------------------#
#
#
#
#
# #------ Second plot: the comparison. ---------------------------------------#
# par(mar=c(4.1,4.6,4.1,2.1))
# plot.new()
# plot.window(xlim=range(montmont),ylim=ylimit,log=plog)
# if (plotgrid){
# abline(v=mplot$levels,h=axTicks(side=2),col=grid.colour,lty="solid")
# }#end if
# if (plotsd){
# polygon(x=mod.x.poly,y=mod.y.poly,col=errcolours[2],angle=angle[2]
# ,density=dens[1],lty="solid",lwd=shwd[1])
# }#end if
# points(x=montmont,y=var.year,col=lcolours[1],lwd=llwd[1],type=ltype
# ,pch=16,cex=1.0)
# points(x=montmont,y=thismean,col=lcolours[2],lwd=llwd[2],type=ltype
# ,pch=16,cex=1.0)
# axis(side=1,at=mplot$levels,labels=mplot$labels,padj=mplot$padj)
# axis(side=2,las=1)
# title(main=letitre,xlab="Time",ylab=ley,cex.main=cex.main)
# box()
# #---------------------------------------------------------------------------#
#
#
#
#
# #----- Close plotting window. ----------------------------------------------#
# if (outform[o] %in% c("x11","quartz")){
# locator(n=1)
# dev.off()
# }else{
# dev.off()
# }#end if
# dummy=clean.tmp()
# #---------------------------------------------------------------------------#
# }#end for outform
# #------------------------------------------------------------------------------#
# }#end for years
# #---------------------------------------------------------------------------------#
# }#end if plotit
# #------------------------------------------------------------------------------------#
# }#end for ncompare
# #---------------------------------------------------------------------------------------#
#
#
#
#
# #---------------------------------------------------------------------------------------#
# # Time series by LU. (tslu folder) #
# #---------------------------------------------------------------------------------------#
# for (v in sequence(ntslu)){
# thistslu = tslu[[v]]
# vnam = thistslu$vnam
# description = thistslu$desc
# unit = thistslu$unit
# plog = thistslu$plog
# plotit = thistslu$plt
#
# #----- Check whether the user wants to have this variable plotted. ------------------#
# if (plotit && any(sellu)){
#
# #---------------------------------------------------------------------------------#
# # Check whether the time series directory exists. If not, create it. #
# #---------------------------------------------------------------------------------#
# outdir = file.path(outpref,"tslu")
# if (! file.exists(outdir)) dir.create(outdir)
# cat(" +",description,"time series for all LUs...","\n")
#
#
#
# #----- Load variable -------------------------------------------------------------#
# thisvar = lu[[vnam]]
# if (plog){
# #----- Eliminate non-positive values in case it is a log plot. ----------------#
# thisvar[thisvar <= 0] = NA
# }#end if
# #---------------------------------------------------------------------------------#
#
# #----- Loop over output formats. -------------------------------------------------#
# for (o in sequence(nout)){
# fichier = file.path(outdir,paste0(vnam,"-",suffix,".",outform[o]))
# if (outform[o] %in% "x11"){
# X11(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "quartz"){
# quartz(width=size$width,height=size$height,pointsize=ptsz)
# }else if(outform[o] %in% "png"){
# png(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent")
# }else if(outform[o] %in% "tif"){
# tiff(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent",compression="lzw")
# }else if(outform[o] %in% "eps"){
# postscript(file=fichier,width=size$width,height=size$height
# ,pointsize=ptsz,paper=size$paper)
# }else if(outform[o] %in% "pdf"){
# pdf(file=fichier,onefile=FALSE
# ,width=size$width,height=size$height,pointsize=ptsz,paper=size$paper)
# }#end if
#
#
# #------------------------------------------------------------------------------#
# # Find the limit, make some room for the legend, and in case the field is #
# # a constant, nudge the limits so the plot command will not complain. #
# #------------------------------------------------------------------------------#
# xlimit = pretty.xylim(u = datum$toyear ,fracexp=0.0,is.log=FALSE)
# ylimit = pretty.xylim(u=thisvar[,sellu],fracexp=0.0,is.log=plog)
# if (plog){
# xylog = "y"
# ydrought = c( exp(sqrt(ylimit[1]^3/ylimit[2]))
# , exp(sqrt(ylimit[2]^3/ylimit[1]))
# )#end c
# }else{
# xylog = ""
# ydrought = c( ylimit[1] - 0.5 * diff(ylimit),ylimit[2] + 0.5 * diff(ylimit) )
# }#end if
# #------------------------------------------------------------------------------#
#
#
#
# #----- Plot settings. ---------------------------------------------------------#
# letitre = paste(description,lieu,sep=" - ")
# ley = desc.unit(desc=description,unit=unit)
# cols = lucols[sellu]
# legs = lunames[sellu]
# #------------------------------------------------------------------------------#
#
#
# #------------------------------------------------------------------------------#
# # Split the plot into two windows. #
# #------------------------------------------------------------------------------#
# par(par.user)
# layout(mat=rbind(2,1),heights=c(5,1))
# #------------------------------------------------------------------------------#
#
#
#
# #------------------------------------------------------------------------------#
# # First plot: legend. #
# #------------------------------------------------------------------------------#
# par(mar=c(0.1,4.6,0.1,2.1))
# plot.new()
# plot.window(xlim=c(0,1),ylim=c(0,1))
# legend( x = "bottom"
# , inset = 0.0
# , legend = legs
# , col = cols
# , lwd = lwidth
# , ncol = min(3,pretty.box(n.sellu)$ncol)
# , title = expression(bold("Land use type"))
# , xpd = TRUE
# , bty = "n"
# )#end legend
# #------------------------------------------------------------------------------#
#
#
#
# #------------------------------------------------------------------------------#
# # Main plot. #
# #------------------------------------------------------------------------------#
# par(mar=c(4.1,4.6,4.1,2.1))
# plot.new()
# plot.window(xlim=xlimit,ylim=ylimit,log=xylog)
# axis(side=1)
# axis(side=2,las=1)
# box()
# title(main=letitre,xlab="Year",ylab=ley,cex.main=0.7)
# if (drought.mark){
# for (n in sequence(ndrought)){
# rect(xleft = drought[[n]][1],ybottom = ydrought[1]
# ,xright = drought[[n]][2],ytop = ydrought[2]
# ,col = grid.colour,border=NA)
# }#end for
# }#end if
# #----- Plot grid. -------------------------------------------------------------#
# if (plotgrid){
# abline(v=axTicks(side=1),h=axTicks(side=2),col=grid.colour,lty="solid")
# }#end if
# #----- Plot lines. ------------------------------------------------------------#
# for (n in sequence(nlu+1)){
# if (sellu[n]){
# lines(datum$toyear,thisvar[,n],type="l",col=lucols[n],lwd=lwidth)
# }#end if
# }#end for
# #------------------------------------------------------------------------------#
#
#
# #----- Close the device. ------------------------------------------------------#
# if (outform[o] %in% c("x11","quartz")){
# locator(n=1)
# dev.off()
# }else{
# dev.off()
# }#end if
# dummy=clean.tmp()
# #------------------------------------------------------------------------------#
# }#end for outform
# #---------------------------------------------------------------------------------#
# }#end if (tseragbpft)
# #------------------------------------------------------------------------------------#
# }#end for tseries
# #---------------------------------------------------------------------------------------#
#
#
#
#
#
# #---------------------------------------------------------------------------------------#
# # Plot disturbance rate by disturbance transition. (distrurb-gyf.pdf file) #
# #---------------------------------------------------------------------------------------#
# if (tserdist && any(seldist)){
# cat(" + Disturbance rate time series for all disturbances...","\n")
# for (o in sequence(nout)){
# fichier = file.path(outpref,paste0("disturb-",suffix,".",outform[o]))
# if (outform[o] %in% "x11"){
# X11(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "quartz"){
# quartz(width=size$width,height=size$height,pointsize=ptsz)
# }else if(outform[o] %in% "png"){
# png(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent")
# }else if(outform[o] %in% "tif"){
# tiff(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent",compression="lzw")
# }else if(outform[o] %in% "eps"){
# postscript(file=fichier,width=size$width,height=size$height
# ,pointsize=ptsz,paper=size$paper)
# }else if(outform[o] %in% "pdf"){
# pdf(file=fichier,onefile=FALSE
# ,width=size$width,height=size$height,pointsize=ptsz,paper=size$paper)
# }#end if
#
# #---------------------------------------------------------------------------------#
# # Find the limit, make some room for the legend, and in case the field is a #
# # constant, nudge the limits so the plot command will not complain. #
# #---------------------------------------------------------------------------------#
# xlimit = pretty.xylim(u=datum$toyear,fracexp=0.0,is.log=FALSE)
# ylimit = NULL
# n = 0
# mylucols = NULL
# mylulegs = NULL
# for (jlu in sequence(nlu)){
# for (ilu in sequence(nlu)){
# n = n + 1
# if (seldist[ilu,jlu]){
# ylimit = c(ylimit,lu$dist[,ilu,jlu])
# mylucols = c(mylucols,distcols [n])
# mylulegs = c(mylulegs,distnames[n])
# }#end if
# }#end for
# }#end for
# ylimit = pretty.xylim(u=ylimit,fracexp=0.0,is.log=FALSE)
# ydrought = c(ylimit[1] - 0.5 * diff(ylimit), ylimit[2] + 0.5 * diff(ylimit))
# #---------------------------------------------------------------------------------#
#
#
#
# #----- Plot settings. ------------------------------------------------------------#
# letitre = paste("Disturbance rates",lieu,sep=" - ")
# #---------------------------------------------------------------------------------#
#
#
# #---------------------------------------------------------------------------------#
# # Split the plot into two windows. #
# #---------------------------------------------------------------------------------#
# par(par.user)
# layout(mat=rbind(2,1),heights=c(5,1))
# #---------------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------------#
# # First plot: legend. #
# #---------------------------------------------------------------------------------#
# par(mar=c(0.1,4.6,0.1,2.1))
# plot.new()
# plot.window(xlim=c(0,1),ylim=c(0,1))
# legend( x = "bottom"
# , inset = 0.0
# , bg = background
# , legend = mylulegs
# , col = mylucols
# , lwd = lwidth
# , ncol = min(3,pretty.box(n)$ncol)
# , title = expression(bold("Transition"))
# , xpd = TRUE
# , bty = "n"
# )#end legend
# #---------------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------------#
# # Main plot. #
# #---------------------------------------------------------------------------------#
# par(mar=c(4.1,4.6,4.1,2.1))
# plot.new()
# plot.window(xlim=xlimit,ylim=ylimit,log=xylog)
# axis(side=1)
# axis(side=2,las=1)
# box()
# title( main = letitre
# , xlab = "Year"
# , ylab = desc.unit(desc="Disturbance rate",unit=untab$oneoyr)
# , cex.main = 0.7
# )#end title
# if (drought.mark){
# for (n in sequence(ndrought)){
# rect(xleft = drought[[n]][1],ybottom = ydrought[1]
# ,xright = drought[[n]][2],ytop = ydrought[2]
# ,col = grid.colour,border=NA)
# }#end for
# }#end if
# #----- Plot grid. ----------------------------------------------------------------#
# if (plotgrid){
# abline(v=axTicks(side=1),h=axTicks(side=2),col=grid.colour,lty="solid")
# }#end if
# #----- Plot lines. ---------------------------------------------------------------#
# n = 0
# for (jlu in sequence(nlu)){
# for (ilu in sequence(nlu)){
# n = n + 1
# if (seldist[ilu,jlu]){
# lines(datum$toyear,lu$dist[,ilu,jlu],type="l"
# ,col=distcols[n],lwd=lwidth)
# }#end if
# }#end for
# }#end for
# #---------------------------------------------------------------------------------#
#
#
# #----- Close the device. ---------------------------------------------------------#
# if (outform[o] %in% c("x11","quartz")){
# locator(n=1)
# dev.off()
# }else{
# dev.off()
# }#end if
# dummy=clean.tmp()
# #---------------------------------------------------------------------------------#
# } #end for outform
# #------------------------------------------------------------------------------------#
# }#end if
# #---------------------------------------------------------------------------------------#
#
#
#
#
# #---------------------------------------------------------------------------------------#
# # Plot the time series diagrams showing annual means. (theme_ymean folder) #
# #---------------------------------------------------------------------------------------#
# cat(" * Plot time series of groups of variables...","\n")
# for (hh in sequence(ntheme)){
#
# #----- Retrieve variable information from the list. ---------------------------------#
# themenow = theme[[hh]]
# vnames = themenow$vnam
# description = themenow$desc
# lcolours = themenow$colour
# llwd = themenow$lwd
# ltype = themenow$type
# plog = themenow$plog
# prefix = themenow$prefix
# group = themenow$title
# unit = themenow$unit
# legpos = themenow$legpos
# plotit = themenow$ymean
# ylimit.fix = themenow$ymean.lim
#
# if (plotit){
#
# #---------------------------------------------------------------------------------#
# # Check whether the time series directory exists. If not, create it. #
# #---------------------------------------------------------------------------------#
# outdir = file.path(outpref,"theme_ymean")
# if (! file.exists(outdir)) dir.create(outdir)
# cat(" +",group,"time series...","\n")
#
#
# #----- Define the number of layers. ----------------------------------------------#
# nlayers = length(vnames)
# #---------------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------------#
# # Find the limit, make some room for the legend, and in case the field is a #
# # constant, nudge the limits so the plot command will not complain. #
# #---------------------------------------------------------------------------------#
# xlimit = pretty.xylim(u=datum$toyear,fracexp=0.0,is.log=FALSE)
# if (any(! is.finite(ylimit.fix))){
# ylimit = NULL
# for (l in sequence(nlayers)){
# thisvar = ymean[[vnames[l]]]
# ylimit = range(c(ylimit,thisvar),na.rm=TRUE)
# }#end for
# ylimit = pretty.xylim(u=ylimit,fracexp=0.0,is.log=plog)
# }else{
# ylimit = ylimit.fix
# }#end if
# if (plog) {
# xylog = "y"
# ydrought = c( exp(sqrt(ylimit[1]^3/ylimit[2]))
# , exp(sqrt(ylimit[2]^3/ylimit[1]))
# )#end c
# }else{
# xylog = ""
# ydrought = c(ylimit[1]-0.5*diff(ylimit),ylimit[2]+0.5*diff(ylimit))
# }#end if
# #---------------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------------#
# # Check if the directory exists. If not, create it. #
# #---------------------------------------------------------------------------------#
#
# #----- Loop over formats. --------------------------------------------------------#
# for (o in sequence(nout)){
# #------ Open file. ------------------------------------------------------------#
# fichier = file.path(outdir,paste0(prefix,"-",suffix,".",outform[o]))
# if (outform[o] %in% "x11"){
# X11(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "quartz"){
# quartz(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "png"){
# png(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent")
# }else if (outform[o] %in% "tiff"){
# tiff(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent",compression="lzw")
# }else if (outform[o] %in% "eps"){
# postscript(file=fichier,width=size$width,height=size$height
# ,pointsize=ptsz,paper=size$paper)
# }else if (outform[o] %in% "pdf"){
# pdf(file=fichier,onefile=FALSE
# ,width=size$width,height=size$height,pointsize=ptsz,paper=size$paper)
# }#end if
# #------------------------------------------------------------------------------#
#
#
#
# #----- Plot settings. ---------------------------------------------------------#
# letitre = paste0(" Time series: ",group,"\n",lieu)
# ley = desc.unit(desc=group,unit=unit)
# #------------------------------------------------------------------------------#
#
#
# #------------------------------------------------------------------------------#
# # Split the plot into two windows. #
# #------------------------------------------------------------------------------#
# par(par.user)
# layout(mat=rbind(2,1),heights=c(5,1))
# #------------------------------------------------------------------------------#
#
#
#
# #------------------------------------------------------------------------------#
# # First plot: legend. #
# #------------------------------------------------------------------------------#
# par(mar=c(0.1,4.6,0.1,2.1))
# plot.new()
# plot.window(xlim=c(0,1),ylim=c(0,1))
# legend( x = "bottom"
# , inset = 0.0
# , legend = description
# , col = lcolours
# , lwd = llwd
# , ncol = min(3,pretty.box(nlayers)$ncol)
# , xpd = TRUE
# , bty = "n"
# )#end legend
# #------------------------------------------------------------------------------#
#
#
#
# #------------------------------------------------------------------------------#
# # Main plot. #
# #------------------------------------------------------------------------------#
# par(mar=c(4.1,4.6,4.1,2.1))
# plot.new()
# plot.window(xlim=xlimit,ylim=ylimit,log=xylog)
# axis(side=1)
# axis(side=2,las=1)
# box()
# title(main=letitre,xlab="Year",ylab=ley,cex.main=0.7)
# if (drought.mark){
# for (n in sequence(ndrought)){
# rect(xleft = drought[[n]][1],ybottom = ydrought[1]
# ,xright = drought[[n]][2],ytop = ydrought[2]
# ,col = grid.colour,border=NA)
# }#end for
# }#end if
# #----- Plot grid. -------------------------------------------------------------#
# if (plotgrid){
# abline(v=axTicks(side=1),h=axTicks(side=2),col=grid.colour,lty="solid")
# }#end if
# #----- Plot lines. ------------------------------------------------------------#
# for (l in sequence(nlayers)){
# thisvar = ymean[[vnames[l]]]
# points(x=datum$toyear,y=thisvar,col=lcolours[l],lwd=llwd[l],type=ltype
# ,pch=16,cex=0.8)
# }#end for
# #------------------------------------------------------------------------------#
#
#
# #----- Close the device. ------------------------------------------------------#
# if (outform[o] %in% c("x11","quartz")){
# locator(n=1)
# dev.off()
# }else{
# dev.off()
# }#end if
# dummy=clean.tmp()
# #------------------------------------------------------------------------------#
# } #end for outform
# #---------------------------------------------------------------------------------#
# }#end if plotit
# #------------------------------------------------------------------------------------#
# }#end for ntser
# #---------------------------------------------------------------------------------------#
#
#
#
#
# #---------------------------------------------------------------------------------------#
# # Plot the climatology of the mean diurnal cycle. (theme_qmean folder) #
# #---------------------------------------------------------------------------------------#
# cat(" * Plot mean diel for groups of variables...","\n")
# for (hh in sequence(ntheme)){
#
# #----- Retrieve variable information from the list. ---------------------------------#
# themenow = theme[[hh]]
# vnames = themenow$vnam
# description = themenow$desc
# lcolours = themenow$colour
# llwd = themenow$lwd
# ltype = themenow$type
# plog = themenow$plog
# prefix = themenow$prefix
# group = themenow$title
# unit = themenow$unit
# legpos = themenow$legpos
# ylimit.fix = themenow$qmean.lim
# plotit = themenow$qmean && plot.ycomp
# if (plog){
# xylog = "y"
# }else{
# xylog = ""
# }#end if
#
# if (plotit){
#
# #---------------------------------------------------------------------------------#
# # Check whether the time series directory exists. If not, create it. #
# #---------------------------------------------------------------------------------#
# outdir = file.path(outpref,"theme_qmean")
# if (! file.exists(outdir)) dir.create(outdir)
# outtheme = file.path(outdir,prefix)
# if (! file.exists(outtheme)) dir.create(outtheme)
# cat(" +",group," diurnal cycle...","\n")
#
#
# #----- Define the number of layers. ----------------------------------------------#
# nlayers = length(vnames)
# xlimit = range(thisday)
# if (any(! is.finite(ylimit.fix))){
# ylimit = NULL
# for (l in sequence(nlayers)) ylimit = c(ylimit,umean[[vnames[l]]])
# ylimit = pretty.xylim(u=ylimit,fracexp=0.0,is.log=plog)
# }else{
# ylimit = ylimit.fix
# }#end if
# #---------------------------------------------------------------------------------#
#
#
# #---------------------------------------------------------------------------------#
# # Loop over all months. #
# #---------------------------------------------------------------------------------#
# yplot = as.numeric(dimnames(umean[[vnames[1]]])[[1]])
# for (yy in seq_along(yplot)){
# cyear = sprintf("%4.4i",yplot[yy])
#
# #------------------------------------------------------------------------------#
# # Check if the directory exists. If not, create it. #
# #------------------------------------------------------------------------------#
#
# #----- Loop over formats. -----------------------------------------------------#
# for (o in sequence(nout)){
# #------ Open file. ---------------------------------------------------------#
# fichier = file.path( outtheme
# , paste0(prefix,"-",cyear,"-",suffix,".",outform[o])
# )#end file.path
# if (outform[o] %in% "x11"){
# X11(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "quartz"){
# quartz(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "png"){
# png(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent")
# }else if (outform[o] %in% "tif"){
# tiff(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent",compression="lzw")
# }else if (outform[o] %in% "eps"){
# postscript(file=fichier,width=size$width,height=size$height
# ,pointsize=ptsz,paper=size$paper)
# }else if (outform[o] %in% "pdf"){
# pdf(file=fichier,onefile=FALSE
# ,width=size$width,height=size$height,pointsize=ptsz,paper=size$paper)
# }#end if
# #---------------------------------------------------------------------------#
#
#
#
# #----- Plot settings. ------------------------------------------------------#
# letitre = paste0(group," - ",lieu,"\n","Mean diurnal cycle - ",cyear)
# ley = desc.unit(desc=group,unit=unit)
# #---------------------------------------------------------------------------#
#
#
# #---------------------------------------------------------------------------#
# # Split the plot into two windows. #
# #---------------------------------------------------------------------------#
# par(par.user)
# layout(mat=rbind(2,1),heights=c(5,1))
# #------------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------#
# # First plot: legend. #
# #---------------------------------------------------------------------------#
# par(mar=c(0.1,4.6,0.1,2.1))
# plot.new()
# plot.window(xlim=c(0,1),ylim=c(0,1))
# legend( x = "bottom"
# , inset = 0.0
# , legend = description
# , col = lcolours
# , lwd = llwd
# , ncol = min(3,pretty.box(nlayers)$ncol)
# , xpd = TRUE
# , bty = "n"
# )#end legend
# #---------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------#
# # Main plot. #
# #---------------------------------------------------------------------------#
# par(mar=c(4.1,4.6,4.1,2.1))
# plot.new()
# plot.window(xlim=xlimit,ylim=ylimit,log=xylog)
# axis(side=1,at=uplot$levels,labels=uplot$labels,padj=uplot$padj)
# axis(side=2,las=1)
# box()
# title(main=letitre,xlab="Year",ylab=ley,cex.main=0.7)
# #----- Plot grid. ----------------------------------------------------------#
# if (plotgrid){
# abline(v=uplot$levels,h=axTicks(side=2),col=grid.colour,lty="solid")
# }#end if
# #----- Plot lines. ---------------------------------------------------------#
# for (l in sequence(nlayers)){
# thisvar = umean[[vnames[l]]]
# thisvar = cbind(thisvar[,ndcycle],thisvar)
# points(x=thisday,y=thisvar[yy,],col=lcolours[l]
# ,lwd=llwd[l],type=ltype,pch=16)
# }#end for
# #---------------------------------------------------------------------------#
#
#
# #----- Close the device. ---------------------------------------------------#
# if (outform[o] %in% c("x11","quartz")){
# locator(n=1)
# dev.off()
# }else{
# dev.off()
# }#end if
# dummy=clean.tmp()
# #---------------------------------------------------------------------------#
# } #end for outform
# #------------------------------------------------------------------------------#
# }#end for pmon
# #---------------------------------------------------------------------------------#
# }#end if plotit
# #------------------------------------------------------------------------------------#
# }#end for ntser
# #---------------------------------------------------------------------------------------#
#
#
#
#
#
# #---------------------------------------------------------------------------------------#
# # Plot the climatology of the soil properties. (soil_ymean folder) #
# #---------------------------------------------------------------------------------------#
# for (v in sequence(nsoilplot)){
#
# #----- Retrieve variable information from the list. ---------------------------------#
# thisclim = soilplot[[v]]
# vnam = thisclim$vnam
# description = thisclim$desc
# unit = thisclim$unit
# vcscheme = thisclim$csch
# pnlog = thisclim$pnlog
# plotit = thisclim$ymean
#
# if (plotit){
#
# #---------------------------------------------------------------------------------#
# # Check if the directory exists. If not, create it. #
# #---------------------------------------------------------------------------------#
# outdir = file.path(outpref,"soil_ymean")
# if (! file.exists(outdir)) dir.create(outdir)
# cat(" + Climatology profile of ",description,"...","\n")
#
# #----- Find the number of rows and columns, and the axes. ------------------------#
# monaxis = sort(unique(datum$year))
# soilaxis = slz
# nmon = length(monaxis)
# nsoil = nzg
#
# #----- Convert the vector data into an array. ------------------------------------#
# vararr = ymean[[vnam]]
#
# #----- Copy the first and the last year to make the edges buffered. --------------#
# first = vararr[1,]
# first = c(first,first[nzg],first[nzg])
#
# last = vararr[nyears,]
# last = c(last[1],last[1],last)
# #---------------------------------------------------------------------------------#
#
#
#
# #----- Bind first and last year to the array, to make the edges buffered. --------#
# varbuff = cbind(vararr[,1],vararr,vararr[,nzg])
# varbuff = rbind(last,varbuff,first)
# #---------------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------------#
# # Expand the month and year axes. #
# #---------------------------------------------------------------------------------#
# yearaxis = c(min(datum$toyear)-1,datum$toyear,max(datum$toyear)+1)
# soilaxis = -log(-1.0 * c( slz[1]*(slz[1]/slz[2])
# , soilaxis
# , slz[nzg]*(slz[nzg]/slz[nzg-1]) ))
#
# if (pnlog){
# vrange = range(varbuff,na.rm=TRUE)
# vlevels = pretty.log(x=vrange,n=ncolshov)
# vnlev = length(vlevels)
# }else{
# vrange = range(varbuff,na.rm=TRUE)
# vlevels = pretty(x=vrange,n=ncolshov)
# vnlev = length(vlevels)
# }#end if
#
# #----- Loop over formats. --------------------------------------------------------#
# for (o in sequence(nout)){
# fichier = file.path(outdir,paste0(vnam,"-",suffix,".",outform[o]))
# if (outform[o] %in% "x11"){
# X11(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "quartz"){
# quartz(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "png"){
# png(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent")
# }else if (outform[o] %in% "tif"){
# tiff(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth,bg="transparent",compression="lzw")
# }else if (outform[o] %in% "eps"){
# postscript(file=fichier,width=size$width,height=size$height
# ,pointsize=ptsz,paper=size$paper)
# }else if (outform[o] %in% "pdf"){
# pdf(file=fichier,onefile=FALSE
# ,width=size$width,height=size$height,pointsize=ptsz,paper=size$paper)
# }#end if
#
# letitre = paste0(description,"\n",lieu)
# ley = desc.unit(desc="Soil depth",unit=untab$m)
# lacle = desc.unit(desc=NULL,unit=unit)
# par(par.user)
# sombreado(x=yearaxis,y=soilaxis,z=varbuff,levels=vlevels,nlevels=vnlev
# ,color.palette=get(vcscheme)
# ,plot.title=title(main=letitre,xlab="Month",ylab=ley,cex.main=0.7)
# ,key.title=title(main=lacle,cex.main=0.8)
# ,key.log=pnlog
# ,plot.axes={axis(side=1)
# axis(side=2,at=zat,labels=znice)
# if (hovgrid){
# abline(h=zat,v=axTicks(1),col=grid.colour,lty="dotted")
# }#end if hovgrid
# }#end plot.axes
# )
#
# if (outform[o] %in% c("x11","quartz")){
# locator(n=1)
# dev.off()
# }else{
# dev.off()
# }#end if
# dummy = clean.tmp()
# } #end for outform
# }#end if plotit
# }#end for nhov
# #---------------------------------------------------------------------------------------#
#
#
#
#
#
#
# #---------------------------------------------------------------------------------------#
# # Bar plot by DBH class. (barplot_dbh folder) #
# #---------------------------------------------------------------------------------------#
# cat(" + Bar plot by DBH classes...","\n")
# pftuse = which(apply(X=szpft$nplant,MARGIN=3,FUN=sum,na.rm=TRUE) > 0.)
# pftuse = pftuse[pftuse != (npft+1)]
# npftuse = length(pftuse)
# pftname.use = pft$name [pftuse]
# pftcol.use = pft$colour[pftuse]
# my.avg = 0
# for (v in sequence(ntspftdbh)){
# #----- Load settings for this variable.----------------------------------------------#
# thisbar = tspftdbh[[v]]
# vnam = thisbar$vnam
# description = thisbar$desc
# unit = thisbar$e.unit
# stacked = thisbar$stack
# plotit = thisbar$bar.plot # && plot.ycomp
# plog = thisbar$plog
# if (plog){
# stacked = FALSE
# xylog = "y"
# }else{
# xylog = ""
# }#end if
# #------------------------------------------------------------------------------------#
#
#
# #------------------------------------------------------------------------------------#
# # Check whether to plot this
# #------------------------------------------------------------------------------------#
# if (plotit){
# cat(" - ",description,"...","\n")
#
#
# #---------------------------------------------------------------------------------#
# # Retrieve the variable, and keep only the part that is usable. #
# #---------------------------------------------------------------------------------#
# thisvnam = szpft[[vnam]]
# thisvnam = thisvnam [,,pftuse]
# thisvnam = thisvnam [,-(ndbh+1),]
#
# thisvnam[is.na(thisvnam)] = 0.
# #---------------------------------------------------------------------------------#
#
#
# #---------------------------------------------------------------------------------#
# # Find the limits for the plots. We use the same axis so it is easier to #
# # compare different times. #
# #---------------------------------------------------------------------------------#
# if (stacked){
# ylimit = c(0,max(apply(X=thisvnam,MARGIN=c(1,2),FUN=sum,na.rm=TRUE)))
# }else{
# ylimit = range(x=thisvnam,na.rm=TRUE)
# }#end if
# ylimit = pretty.xylim(u=ylimit,fracexp=0.0,is.log=plog)
# #---------------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------------#
# # Check if the directory exists. If not, create it. #
# #---------------------------------------------------------------------------------#
# barplotdir = file.path(outpref,"barplot_dbh")
# if (! file.exists(barplotdir)) dir.create(barplotdir)
# outdir = file.path(barplotdir,vnam)
# if (! file.exists(outdir)) dir.create(outdir)
# #---------------------------------------------------------------------------------#
#
#
#
# #---------------------------------------------------------------------------------#
# # Loop over all possible months. #
# #---------------------------------------------------------------------------------#
# for (y in sequence(nyears)){
#
# #----- Find which year we are plotting. ---------------------------------------#
# cyear = sprintf("%4.4i",datum$toyear[y])
# yy = as.numeric(cyear)
# #------------------------------------------------------------------------------#
#
#
# #----- Loop over output formats. ----------------------------------------------#
# for (o in sequence(nout)){
# #------ Open the plot. -----------------------------------------------------#
# fichier = file.path( outdir
# , paste0(vnam,"-",cyear,"-",suffix,".",outform[o])
# )#end file.path
# if (outform[o] %in% "x11"){
# X11(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "quartz"){
# quartz(width=size$width,height=size$height,pointsize=ptsz)
# }else if (outform[o] %in% "png"){
# png(filename=fichier,width=size$width*depth,height=size$height*depth
# ,pointsize=ptsz,res=depth)
# }else if (outform[o] %in% "eps"){
# postscript(file=fichier,width=size$width,height=size$height
# ,pointsize=ptsz,paper=size$paper)
# }else if (outform[o] %in% "pdf"){
# pdf(file=fichier,onefile=FALSE
# ,width=size$width,height=size$height,pointsize=ptsz,paper=size$paper)
# }#end if
# #---------------------------------------------------------------------------#
#
# #----- Split plotting area into two. ---------------------------------------#
# par(par.user)
# layout(mat=rbind(2,1),heights=c(5,1))
# #---------------------------------------------------------------------------#
#
#
#
# #----- Plot the legend. ----------------------------------------------------#
# par(mar=c(0.1,4.6,0.1,2.1))
# plot.new()
# plot.window(xlim=c(0,1),ylim=c(0,1))
# legend( x = "bottom"
# , inset = 0.0
# , legend = pftname.use
# , fill = pftcol.use
# , ncol = min(3,pretty.box(n.selpft)$ncol)
# , title = expression(bold("Plant functional type"))
# , cex = cex.ptsz
# , bg = background
# , xpd = TRUE
# )
# #---------------------------------------------------------------------------#
#
#
# #------ Set up the title and axis labels. ----------------------------------#
# letitre = paste0(lieu,"\n",description," - Year : ",cyear)
# lexlab = "DBH Classes"
# leylab = desc.unit(desc=description,unit=unit)
# #---------------------------------------------------------------------------#
#
#
# #------ Correct data matrix for log plot (0s create problems, NAs not) -----#
# if (xylog == "y"){
# temp = thisvnam[y,,]
# temp[temp==0] = NA
# thisvnam[y,,] = temp
# }
# #---------------------------------------------------------------------------#
#
#
# #------- Plot all monthly means together. ------------------------------------#
# par(mar=c(4.1,4.6,4.1,2.1))
# barplot(height=t(thisvnam[y,,]),names.arg=dbhnames[1:ndbh],width=1.0,
# main=letitre,xlab=lexlab,ylab=leylab,ylim=ylimit * 1.05,
# legend.text=FALSE, beside=(! stacked),col=pftcol.use,log=xylog,
# border=grey.fg,xpd=FALSE,cex.main=cex.main,las=1)
# if (plotgrid & (! stacked)){
# xgrid=0.5+(1:ndbh)*(1+npftuse)
# abline(v=xgrid,col=grid.colour,lty="solid")
# }
# box()
# #---------------------------------------------------------------------------#
# if( v == 24 & plot.nplant.hystogram & y >= (nyears-100)){
# # dev.off()
# xylog = ""
# ylimit=c(0,140.0)
#
# mytable = thisvnam[y,,]
# mytable[is.na(mytable)] = 0.0
# my.avg=my.avg + mytable
# if (y == nyears){
#
# #multiply by ha conversion and divide by 100 (=nyears) (hence *100)
# my.avg = my.avg*100
#
# exp_lianas_dist = c(1000,117.1,53.594,28.79,15.454,
# 12.631,9.331,5.552,5.357,3.728,
# 2.816,1.188,0.754,0.573,0.5,0.5
# ,0.5)
# exp_lianas_dist_err = exp_lianas_dist * 0.3 + 2
# exp_lianas_undist = c(1000,52,31.831,20.801,10.735,
# 10.29,5.996,4.11,2.701,2.258,
# 1.814,1.369,0.5,0.5,0.5,0.5
# ,0.5)
# exp_lianas_undist_err = exp_lianas_undist * 0.4 + 2
#
# my.avg[,1] = exp_lianas_dist
# my.avg[,2] = exp_lianas_undist
#
# my.avg=my.avg[,-3]
#
# er_mat = matrix(c(exp_lianas_dist_err, exp_lianas_undist_err), nrow = 17)
# er_mat = cbind (er_mat,NA)
# er_mat [2,2] = er_mat[2,2] * 1.3
#
#
# fichier=paste(outroot,"/paracou/yearly/barplot_dbh/nplant/average.pdf",sep='')
# pdf(file=fichier,onefile=FALSE
# ,width=8.8,height=6.8,pointsize=16,paper="special")
# #----- Split plotting area into two. ---------------------------------------#
# par(par.user)
# layout(mat=rbind(2,1),heights=c(5,1))
# #---------------------------------------------------------------------------#
#
#
#
# #----- Plot the legend. ----------------------------------------------------#
# par(mar=c(0.1,4.6,0.1,2.1))
# plot.new()
# plot.window(xlim=c(0,1),ylim=c(0,1))
# legend( x = "bottom"
# , inset = 0.0
# , legend = c("Disturbed plots","Undistrurbed plots","Model")
# , fill = c("#D3D3D3","#838B8B","#1E64C8")
# , ncol = min(3,pretty.box(n.selpft)$ncol)
# , title = expression(bold("Type"))
# , cex = cex.ptsz
# , bg = background
# , xpd = TRUE
# )
# #---------------------------------------------------------------------------#
#
#
# #------ Set up the title and axis labels. ----------------------------------#
# lexlab = "DBH Classes"
# leylab = desc.unit(desc="Plant density",unit="p*l*a*n*t^phantom(1)*h*a^{-1}")
# #---------------------------------------------------------------------------#
#
# pftcol.use = c("#D3D3D3","#838B8B","#1E64C8")
#
# #------- Plot all monthly means together. ------------------------------------#
# par(mar=c(4.1,4.6,4.1,2.1))
# mynames = c("<2", "2-3", "3-4", "4-5", "5-6", "6-7", "7-8",
# "8-9", "9-10", "10-11", "11-12", " 12-13", "13-14",
# "14-15", "15-16", "16-20", ">20")
# lexlab="DBH Classes [cm]"
# barx = barplot(height=t(my.avg),names.arg=mynames,width=1.0,
# main="Diameter distribution",xlab=lexlab,ylab=leylab,ylim=ylimit,
# legend.text=FALSE, beside=(! stacked),col=pftcol.use,log="",
# border=grey.fg,xpd=FALSE,cex.main=cex.main,las=1)
# if (plotgrid & (! stacked)){
# npftuse=3
# xgrid=0.5+(1:ndbh)*(1+npftuse)
# abline(v=xgrid,col=grid.colour,lty="solid")
# }
# arrows(barx, y0=t(my.avg)-t(er_mat)/2 , y1=t(my.avg)+t(er_mat) / 2, barx, length=0)
# # arrows(barx, y1=t(my.avg)+t(er_mat), barx, height=t(my.avg), angle=90, code=1, length=0)
# box()
# #---------------------------------------------------------------------------#
#
# }
# }
#
#
#
# #---------------------------------------------------------------------------#
# # Close the device. #
# #---------------------------------------------------------------------------#
# if (outform[o] %in% c("x11","quartz")){
# locator(n=1)
# dev.off()
# }else{
# dev.off()
# }#end if
# dummy = clean.tmp()
# #---------------------------------------------------------------------------#
# } #end for outform
# #------------------------------------------------------------------------------#
# }#end for
# #---------------------------------------------------------------------------------#
# }#end if
# #------------------------------------------------------------------------------------#
# }#end for
# #---------------------------------------------------------------------------------------#
#
#---------------------------------------------------------------------------------------#
# Plot the 3-D size and age structure of various variables. (sas folder) #
#---------------------------------------------------------------------------------------#
#for (v in sequence(ntspftdbh)){
v=69
#----- Retrieve variable information from the list. ---------------------------------#
thissas = tspftdbh[[v]]
vnam = thissas$vnam
description = thissas$desc
unit = thissas$i.unit
plotit = thissas$sas
plog = thissas$plog
#----- If this variable is to be plotted, then go through this if block. ------------#
if (plotit){
cat(" + Size and age structure plot: ",description,"...","\n")
#---------------------------------------------------------------------------------#
# Check if the directory exists. If not, create it. #
#---------------------------------------------------------------------------------#
sasdir = file.path(outpref,"sas")
if (! file.exists(sasdir)) dir.create(sasdir)
outdir = file.path(sasdir,vnam)
if (! file.exists(outdir)) dir.create(outdir)
#---------------------------------------------------------------------------------#
#----- Load this list into "thislist". -------------------------------------------#
varco = cohort[[vnam]]
#---------------------------------------------------------------------------------#
#---------------------------------------------------------------------------------#
# Loop over all times. #
#---------------------------------------------------------------------------------#
for (ww in names(cohort$age)){
#----- Find which year we are plotting. ---------------------------------------#
cmonth = substring(ww,7,8)
thisyear = substring(ww,2,5)
mm = as.numeric(cmonth)
yy = as.numeric(thisyear)
thisyear = as.character(yy - 350)
#------------------------------------------------------------------------------#
mxww = numeric()
myww = numeric()
mzww = numeric()
mpchww = numeric()
mcolww = numeric()
mcexww = numeric()
#----- Retrieve variable list, age, DBH, and PFT for this year. ---------------#
ageww = cohort$age [[ww]]
if (any(ageww <= 0,na.rm=TRUE)){
minww = min(ageww,na.rm=TRUE)
ageww = ageww - minww + 0.01
}#end if
dbhww = cohort$dbh [[ww]]
pftww = cohort$pft [[ww]]
varww = varco [[ww]]
popww = cohort$nplant[[ww]] * cohort$area[[ww]]
#------------------------------------------------------------------------------#
#------------------------------------------------------------------------------#
# We only plot the SAS figures when the polygon is not an absolute desert. #
#------------------------------------------------------------------------------#
if (any (! is.na(varww))){
#---------------------------------------------------------------------------#
# Find the range. If the user wants the range to be fixed, then use #
# the global range, otherwise, simply use the range for this year. #
#---------------------------------------------------------------------------#
if (sasfixlimits){
xlimit = pretty.xylim(u=unlist(cohort$age),fracexp=0.0,is.log=TRUE )
ylimit = pretty.xylim(u=unlist(cohort$dbh),fracexp=0.0,is.log=FALSE)
zlimit = pretty.xylim(u=unlist(varco) ,fracexp=0.0,is.log=plog )
popmin = min(unlist(cohort$nplant * cohort$area), na.rm=TRUE)
popmax = max(unlist(cohort$nplant * cohort$area), na.rm=TRUE)
}else{
xlimit = pretty.xylim(u=ageww ,fracexp=0.0,is.log=TRUE )
ylimit = pretty.xylim(u=dbhww ,fracexp=0.0,is.log=FALSE)
zlimit = pretty.xylim(u=varww ,fracexp=0.0,is.log=plog )
popmin = min(popww ,na.rm=TRUE)
popmax = max(popww ,na.rm=TRUE)
}#end if
#---------------------------------------------------------------------------#
#----- Define the scale-dependent population size. -------------------------#
cexww = cexmin + (cexmax - cexmin) * log(popww/popmin) / log(popmax/popmin)
#---------------------------------------------------------------------------#
#----- Define the floor location. ------------------------------------------#
if ((zlimit[1] > 0) != (zlimit[2] > 0)){
floor3d = 0.
}else if (zlimit[1] > 0){
floor3d = zlimit[1]
}else{
floor3d = zlimit[2]
}#end if
#---------------------------------------------------------------------------#
xlimit = c(0.1, 300)
ylimit = c(0.01, 130)
zlimit = c(0.4999, 34.999)
xlimit.stat = c(-2.5, 300)
ylimit.stat = c(0.01, 130)
zlimit.stat = c(0.4999, 34.999)
floor3d = 0.0
#----- Define the grid information for the 3-D plot. -----------------------#
xlabels = pretty.log(xlimit,n=5)
ylabels = pretty(ylimit,n=5)
zlabels = if(plog){pretty.log(zlimit,n=5)}else{pretty(zlimit,n=5)}
xat = log(xlabels)
yat = ylabels
zat = if(plog){log(zlabels)}else{zlabels}
xlimit = range(x=xat)
ylimit = range(x=yat)
zlimit = range(x=zat)
xfloor = seq(from=xlimit[1],to=xlimit[2],length.out=16)
yfloor = seq(from=ylimit[1],to=ylimit[2],length.out=16)
zfloor = matrix(floor3d,nrow=length(xfloor),ncol=length(yfloor))
#---------------------------------------------------------------------------#
#----- Expand the lines to make the lollipops. -----------------------------#
ncohnow = length(varww)
ageww = rep(ageww,each=3)
dbhww = rep(dbhww,each=3)
pftww = rep(pftww,each=3)
varww = as.vector( rbind( rep(floor3d,times=ncohnow)
, varco[[ww]]
, rep(NA,times=ncohnow)
)#end rbind
)#end as.vector
xww = log(ageww)
yww = dbhww
zww = if(plog){log(varww)}else{varww}
pchww = rep(c(NA,16,NA),times=ncohnow)
cexww = rep(cexww,each=3)
colww = pft$colour[pftww]
pftin = sort(unique(cohort$pft[[ww]]))
colleg = pft$colour[pftin]
pftleg = pft$name [pftin]
#---------------------------------------------------------------------------#
#---------------------------------------------------------------------------#
# Plot annotation. #
#---------------------------------------------------------------------------#
letitre = paste(description," - ",lieu
,"\n Time :",mlist[mm],"/",thisyear,sep=" ")
lexlab = desc.unit(desc="Gap age",unit=untab$yr)
leylab = desc.unit(desc="DBH",unit=untab$cm)
lezlab = desc.unit(desc=description,unit=unit)
#---------------------------------------------------------------------------#
#----- Loop over output formats. -------------------------------------------#
for (o in sequence(nout)){
#----- Open file. -------------------------------------------------------#
fichier = file.path( outdir
, paste0( vnam,"-",thisyear,"-",cmonth,"-",suffix
,".",outform[o]
)#end paste0
)#end file.path
if (outform[o] %in% "x11"){
X11(width=size$width,height=size$height,pointsize=ptsz)
}else if (outform[o] %in% "quartz"){
quartz(width=size$width,height=size$height,pointsize=ptsz)
}else if(outform[o] %in% "png"){
png(filename=fichier,width=size$width*depth,height=size$height*depth
,pointsize=ptsz,res=depth,bg="transparent")
}else if(outform[o] %in% "tif"){
tiff(filename=fichier,width=size$width*depth,height=size$height*depth
,pointsize=ptsz,res=depth,bg="transparent",compression="lzw")
}else if(outform[o] %in% "eps"){
postscript(file=fichier,width=size$width,height=size$height
,pointsize=ptsz,paper=size$paper)
}else if(outform[o] %in% "pdf"){
pdf(file=fichier,onefile=FALSE,width=size$width,height=size$height
,pointsize=ptsz,paper=size$paper)
}#end if
#------------------------------------------------------------------------#
#----- Split the domain into 2. -----------------------------------------#
par(par.user)
layout(mat=rbind(2,1),heights=c(5,1))
#------------------------------------------------------------------------#
#------------------------------------------------------------------------#
# Plot legend. #
#------------------------------------------------------------------------#
par(mar=c(0.1,0.1,0.1,0.1))
plot.new()
plot.window(xlim=c(0,1),ylim=c(0,1))
legend( x = "center"
, inset = 0.0
, legend = pftleg
, fill = colleg
, ncol = min(4,pretty.box(length(pftleg))$ncol)
, title = expression(bold("Plant functional type"))
, cex = cex.ptsz
, xpd = TRUE
, bty = "n"
)#end legend
#------------------------------------------------------------------------#
#------------------------------------------------------------------------#
# Plot the 3-D plot. #
#------------------------------------------------------------------------#
par(mar=c(1.1,1.1,4.1,1.1))
pout = perspx( x = xfloor
, y = yfloor
, z = zfloor
, xlim = xlimit
, ylim = ylimit
, zlim = zlimit
, theta = theta
, phi = phi
, col = gcol
, expand = expz
, ticktype = "detailed"
, border = NA
, shade = shade
, ltheta = ltheta
, main = letitre
, cex.main = 0.8*cex.ptsz
, axes = FALSE
)#end perspx
#----- Add axes. --------------------------------------------------------#
paxis3d(edge="X--",pmat=pout,at=xat,cex=0.9*cex.ptsz,labels=xlabels)
paxis3d(edge="Y--",pmat=pout,at=yat,cex=0.9*cex.ptsz,labels=ylabels)
paxis3d(edge="Z-+",pmat=pout,at=zat,cex=0.9*cex.ptsz,labels=zlabels)
mtext3d(edge="X--",pmat=pout,labels=lexlab,cex=cex.ptsz,srt=theta+90)
mtext3d(edge="Y--",pmat=pout,labels=leylab,cex=cex.ptsz,srt=theta)
mtext3d(edge="Z-+",pmat=pout,labels=lezlab,cex=cex.ptsz,srt=-75)
#------------------------------------------------------------------------#
for (i in sequence(length(xww))){
if (xww[i] %wr% xlimit.stat && yww[i] %wr% ylimit.stat){
mxww = c(mxww , xww[i])
myww = c(myww , yww[i])
mzww = c(mzww , zww[i])
mpchww = c(mpchww, pchww[i])
mcolww = c(mcolww, colww[i])
mcexww = c(mcexww, cexww[i])
}
}
if(v==69){
#----- Add the cohorts. -------------------------------------------------#
lines (trans3d(x=mxww,y=myww,z=mzww,pmat=pout),type="l",col=grey.fg,lwd=2)
points(trans3d(x=mxww,y=myww,z=mzww,pmat=pout),type="p",pch=mpchww
,col=mcolww,cex=mcexww)
#------------------------------------------------------------------------#
}
#----- Add the cohorts. -------------------------------------------------#
# lines (trans3d(x=xww,y=yww,z=zww,pmat=pout),type="l",col=grey.fg,lwd=2)
# points(trans3d(x=xww,y=yww,z=zww,pmat=pout),type="p",pch=pchww
# ,col=colww,cex=cexww)
#------------------------------------------------------------------------#
#----- Close the device. ------------------------------------------------#
if (outform[o] %in% c("x11","quartz")){
locator(n=1)
dev.off()
}else{
dev.off()
}#end if
dummy = clean.tmp()
#------------------------------------------------------------------------#
}#end for outform
#---------------------------------------------------------------------------#
}#end if is.na(varww)
#------------------------------------------------------------------------------#
}#end for nameco
#---------------------------------------------------------------------------------#
}#end if
#------------------------------------------------------------------------------------#
#}#end for npsas
#---------------------------------------------------------------------------------------#
#}#end for places
#q("no")
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.