########################################################################################
# Spinup analysis script
# see iland.boku.ac.at/spinup
# This is an example script that works together with the spinup.js/spinuplib.js and iLand.
# Werner Rammer, 2017
########################################################################################

library(stringr)
library(dplyr)
library(tidyr)


############# Inits ##############

cols.species <- c("abal"="#088A29",
                  "acca"="#F3F781",
                  "acpl"="#86B404", 
                  "acps"="#58FAD0", 
                  "algl"="#61210B", 
                  "alin"="#A4A4A4", 
                  "alvi"="#0B3B17", 
                  "bepe"="#2E64FE", 
                  "cabe"= "#F7BE81", 
                  "casa"="#A9F5A9", 
                  "coav"="#58D3F7",
                  "fasy"="#2EFE2E", 
                  "frex"="#FF0000", 
                  "lade"="#8A4B08", 
                  "piab"="#FFFF00", 
                  "pice"="#FF8000", 
                  "pini"="#610B21", 
                  "pisy"="#B18904", 
                  "poni"="#000000", 
                  "potr"="#610B5E", 
                  "psme"="#DF0174", 
                  "qupe"="#F6CED8", 
                  "qupu"="#FA5882", 
                  "quro"="#B40431", 
                  "rops"="#F781D8", 
                  "saca"="#F5ECCE", 
                  "soar"="#E6E0F8", 
                  "soau"="#B40404", 
                  "tico"="#9F81F7",
                  "tipl"="#8000FF", 
                  "ulgl"="#DF7401", "rest"="#202020")

### target values are extracted from the abe-csv file (but could be )
targets <- read.csv("D:/Projects/RESIN/Stubaital/Stubaital_Sim/abe/abe.csv")
target.rel <- targets[, c('id', 'piab', 'abal', 'lade', 'pisy', 'pice', 'broadleaf')]
names(target.rel)[7] <- "rest"
head(target.rel)

########### Analyze a spinup run ##############


# load the log-files created by the iLand spinup scripts:
dbhdist <- read.csv("D:/Projects/RESIN/Stubaital/Stubaital_Sim/temp/dbhdist.txt", sep=";")
spinuplog <- read.csv("D:/Projects/RESIN/Stubaital/Stubaital_Sim/temp/spinup.log", sep=";", stringsAsFactors = FALSE)
stands <- read.csv("D:/Projects/RESIN/Stubaital/Stubaital_Sim/temp/stands.txt", sep=";")


# specify name/path of the output PDF
pdf("spinup_1.11_full.pdf")

par(mfrow=c(1,1))
barplot(as.matrix(dbhdist[, 2:22]), col=cols.species[match(dbhdist[,1], names(cols.species))], main="final dbh distribution")
legend(x="topright", legend=dbhdist[,1], ncol=2, fill=cols.species[match(dbhdist[,1], names(cols.species))])

sl.eval <- spinuplog[grepl("^EVAL", spinuplog$text), ]
vol.table <- as.data.frame(cbind(sl.eval$year, sl.eval$stand, t( sapply( str_match_all(sl.eval$text, "([0-9.]+)"), function(x) as.numeric(x[,1])) ) ))
names(vol.table) <- c("year","stand", "plan", "current")

### species shares of the individual stands (based on basal area)
# in "long" format
sl.shares <- spinuplog[grepl("^share:", spinuplog$text),]
s <- list()
for (i in 1:length(sl.shares$year)) {
  df <-  str_match_all(sl.shares$text[i], "([a-z]+): ([0-9.]+)")[[1]][,2:3]
  if (length(df)>0) {
    if (is.matrix(df))
      s[[i]] <- data.frame(stand=sl.shares$stand[i], year=sl.shares$year[i], df)
    else
      s[[i]] <- data.frame(stand=sl.shares$stand[i], year=sl.shares$year[i], t(df))
  }
}
#rm(sshare)
sshare <- do.call(rbind, s)
names(sshare)[3:4] <- c("species", "share")
sshare$species <- factor( sshare$species, levels=c("piab","abal","lade", "pisy","pice", "rest"))
sshare$species[is.na(sshare$species)] <- 'rest'
sshare$share <- as.numeric( levels(sshare$share)[sshare$share] )
sshare <- sshare %>% group_by(stand,year,species) %>% summarise(share=sum(share))
# convert to wide format
sshare <- spread(sshare, key=species, value=share, fill=0 )

sl.save <- spinuplog[grepl("^save:", spinuplog$text), ]
sl.score <- spinuplog[grepl("^score:", spinuplog$text), ]
## fix order of 

par(mfrow=c(3,3))
for(standid in unique(sshare$stand)) {
  tm <- (t(sshare[sshare$stand==standid,3:8])) 
  tm <- cbind(tm, t(target.rel[target.rel$id==standid, match(rownames(tm), names(target.rel ))]))  ## add a column...
  tn <- c(sshare$year[sshare$stand==standid], rep("target", length(target.rel$id[target.rel$id==standid]) ) )
  
  # add ** to those years where the stand was saved
  tn[tn %in% sl.save$year[sl.save$stand==standid] ] <- paste("*", tn[tn %in% sl.save$year[sl.save$stand==standid] ], "*", sep="")
  l <- strsplit(sl.score$text[sl.score$stand==standid], " ")
  barplot( tm, col=cols.species[match(rownames(tm), names(cols.species))], main=standid, names.arg=tn)
  matplot(t(rbind(sapply(l, function(x) as.numeric(x[2])), sapply(l, function(x) as.numeric(x[3])) )), type="l", ylab="score" )
  legend("topleft", legend = c("score", "species"), col=1:2, lwd=1)
  vtarget <- max( targets$basalarea[targets$id==standid] )
  plot( vtarget-sapply(l, function(x) as.numeric(x[4])), type="b", ylim=c(0,70), ylab="volume", main="basal area" )
  abline(a=vtarget, b=0)
}


ages <- targets$age[match(vol.table$stand, targets$id)]
ages <- cut(ages, breaks=c(0,40,80,120,160, 200))
par(mfrow=c(2,2))
boxplot(vol.table$current-vol.table$plan~ ages, main="Delta m2", xlab="age", ylim=c(-20,20), xlab="Delta (sim-plan)")
abline(a=0,b=0)
boxplot(vol.table$current, vol.table$plan, names=c("Simulated", "Plan"), main="Standing volume (m3/ha)")
boxplot(vol.table$current~ ages, main="basal area simulated", xlab="age", ylim=c(0,80))
boxplot(vol.table$plan~ ages, main="basal area planned", xlab="age", ylim=c(0,80))

  
par(mfrow=c(1,1))
plot(vol.table$current ~ vol.table$plan, pch=20, col=rgb(0,0,0,.5))
abline(a=0,b=1)
summary(lm(vol.table$current ~ vol.table$plan))



hist(sl.save$year, main="# of saved stands", xlab="year")
# volume / basalarea / landscape
stands %>% summarise(volume=sum(volume)/sum(area)*10000, basalArea=sum(basalArea)/sum(area)*10000,stems=sum(stems)/sum(area)*10000)

plotDbhDists <- function(path_filter) {
  files <- dir(path_filter, "dbh*")
  all.vol <- data.frame()
  for (f in files) {
    yr = str_match(f, "_([0-9]+)")[2]
    if (is.na(yr))
      yr="Latest"
    dbhdist <- read.csv(paste(path_filter, "/", f, sep=""), sep=";")
    barplot(as.matrix(dbhdist[, 2:22]), col=cols.species[match(dbhdist[,1], names(cols.species))], main=paste("year",yr))
    #legend(x="topright", legend=dbhdist[,1], ncol=2, fill=cols.species[match(dbhdist[,1], names(cols.species))])
    if (dim(all.vol)[1]==0)
      all.vol <- data.frame(species=dbhdist$species)
    all.vol[paste("yr", yr)] <-  dbhdist$volume
  }
  barplot(as.matrix(all.vol[, 2:dim(all.vol)[2]]), col=cols.species[match(all.vol[,1], names(cols.species))], main="species shares", ylab="volume (m3/ha)")
}

par(mfrow=c(2,2))
plotDbhDists("D:/Projects/RESIN/Stubaital/Stubaital_Sim/temp")

plotStandDev <- function(path_filter) {
  files <- dir(path_filter, "stands*")
  all.stands <- data.frame()
  for (f in files) {
    yr = str_match(f, "_([0-9]+)")[2]
    if (is.na(yr))
      yr="Latest"
    stands <- read.csv(paste(path_filter, "/", f, sep=""), sep=";")
    all.stands <- rbind(all.stands, cbind(stands, yr=yr))
    # add "empty" stands
    #all.stands <- rbind(all.stands, data.frame(standId=targets.0vol$stand_id, area=targets.0vol$area*10000, basalArea=0, volume=0, stems=0, yr=yr))
  }
  boxplot(all.stands$basalArea/all.stands$area*10000 ~ all.stands$yr, ylab="basalarea (m2/ha)")
  boxplot(all.stands$volume/all.stands$area*10000 ~ all.stands$yr, ylab="standing volume (m3/ha)")
  all.stands
}

par(mfrow=c(1,1))
all.stands <- plotStandDev("D:/Projects/RESIN/Stubaital/Stubaital_Sim/temp")


sum.stands <- all.stands %>% group_by(yr)  %>% summarise(volume=sum(volume)/sum(area)*10000, basalArea=sum(basalArea)/sum(area)*10000,stems=sum(stems)/sum(area)*10000)

par(mfrow=c(1,1))
barplot(sum.stands$volume, names.arg=sum.stands$yr, main="volume")
barplot(sum.stands$basalArea, names.arg=sum.stands$yr, main="basal area")
barplot(sum.stands$stems, names.arg=sum.stands$yr, main="Stems (n/ha)")


asl <- all.stands[all.stands$yr=="Latest",]
asl <- left_join(asl, targets, by=c("standId"="id")) 
asl <- colSums(asl[,11:16]*(asl$area/10000))
asl <- asl/ sum(asl, na.rm=T)
names(asl)[names(asl) == "broadleaf"] <- "acps"

tm <- full_join(data.frame(species=names(asl), share=asl), dbhdist[,c("species", "basalArea")])
tm$basalArea <- tm$basalArea/sum(tm$basalArea)
barplot(as.matrix(tm[,2:3]), col=cols.species[match(tm$species, names(cols.species))], names.arg=c("target", "simulated"), ylab="species share (basal area)")


dev.off()