<!--
%\VignetteEngine{litedown::vignette}
%\VignetteIndexEntry{Adiele}
%\VignetteEncoding{UTF-8}
-->

## LINTULcassava model

This document shows how to run the LINTUL-Cassava model and how to produce the figures shown in 

Adiele, J.G., A.G.T. Schut, R.P.M. van den Beuken, K.S. Ezui, P. Pypers, A.O. Ano, C.N. Egesi & K.E. Giller (2021). A recalibrated and tested LINTUL-Cassava simulation model provides insight into the high yield potential of cassava under rainfed conditions. European Journal of Agronomy 124: 126242.

### Run model 

We first run the model for three locations and two years

```{r run_model}
library(LINTULcassava)

#parameters defined by Ezui et al. and estimated by Adiele et al. EJA 2021
crop <- LC_crop("Adiele")

#model run combinations
x <- expand.grid(watlim=c(FALSE, TRUE), site=c("Edo", "CRS", "Benue"), 
                 year = c(2016, 2017), stringsAsFactors=FALSE)

mod <- lapply(1:nrow(x), function(i)  {
    p <- Adiele(x$site[i], x$year[i])
    data.frame(Location=x$site[i], 
        LINTCAS(p$weather, crop, p$soil, p$management, 
			control=c(p$control, water_limited=x$watlim[i]))
    )
})
names(mod) <- apply(x, 1, function(i) paste0(c(ifelse(i[1], "wlim", "pot"), i[2:3]), collapse="_"))
```


### Figure 1A

Helper function to read example data
```{r helper}
lcp <- function(f) file.path(system.file(package="LINTULcassava"), f)
```

Read data and fit regression model

```{r anova}

#Data measured with an AccuPar below and above canopy
DATA_ACCUPAR <- read.csv(lcp("ex/AccuPAR_and_caculated_LAI_k.csv"))
Modk <- lm(ln_FLI ~ -1 + LAI, data=DATA_ACCUPAR)
print(summary(Modk))
print(anova(Modk))
```

Make figure

```{r figure1A, fig.width=7, fig.height=7}

Fig_1A <- function(DATA, model){
  par(oma=c(2,2,2,2),mar=c(5,5,0,0))
  #Modk <- lm(ln_FLI ~ -1 + LAI, data=DATA)
  
  uLY <- unique(DATA[,"Loc_yr"])
  i <- DATA$Loc_yr == "Benue_2017"
  plot(DATA[i,"LAI"], DATA[i,"ln_FLI"], cex.lab=1.5, ylab = '-ln (1-fLI)',
       xlab = 'Leaf area index (m²/m²)', xlim=c(0,7), ylim=c(0,5), pch=15)
  i <- DATA$Loc_yr=="Benue_2016"
  points(x=DATA[i,"LAI"], y=DATA[i,"ln_FLI"], pch=0)
  i <- DATA[,"Loc_yr"]=="Edo_2016"
  points(x=DATA[i,"LAI"], y=DATA[i,"ln_FLI"], pch=1)
  i <- DATA[,"Loc_yr"]=="Edo_2017"
  points(x=DATA[i,"LAI"], y=DATA[i,"ln_FLI"], pch=16)

  legend("bottomright", c("Benue, 2016" ,"Benue, 2017","Edo, 2016", "Edo, 2017"),
          pch = c(0, 15, 1, 16), bty="n") 
  abline(model)
  legend('topleft', c('R² = 0.99', 'y = 0.67x'), bty = 'n')
}

Fig_1A(DATA = DATA_ACCUPAR, Modk)
```

### Figure 1

```{r figure1, fig.width=7, fig.height=7}
Fig_1 <- function(DATA, model){

  par(oma=c(2,2,2,2),mar=c(5,5,0,0))

  plot(DATA[,"LAI"], DATA[,"ln_FLI"], cex.lab=1.5,  
       xlab = 'Leaf area index (m²/m²)', ylab = '-ln (1-fLI)', 
       xlim = c(0,7), ylim = c(0,5))

  abline(model)
  legend('topleft', c('R² = 0.99', 'y = 0.67x'), bty = 'n')
}

Fig_1(DATA = DATA_ACCUPAR, Modk)
```

### Figure 2


```{r fig2}

DAP_vs_IPAR <- function(d, loc, year, leg=FALSE) {
    ly <- d$Location==loc & DATA$Year==year 
    S1 <- d[ly & DATA$Treatment=="NfPfK1", ]
    S2 <- d[ly & DATA$Treatment=="NfPfKf", ]
    S3 <- d[ly & DATA$Treatment=="NfPfKfMn", ]
    LOCYR=paste(loc, year, collapse=", ")
    LEG=leg 
    LEGTEXT=c('NfPfK240', 'NfPfKf', 'NfPfKfMN')
  
    red <- rgb(1, 0, 0, alpha = 0.8)
    blue <- rgb(0, 0, 1, alpha = 0.8)
  
  plot(S1$DAP, S1$fPAR, type="l", xlab = "Days after planting", 
        ylab = "IPAR, f",  ylim = c(0,1), xlim=c(0,450), lwd = 1.5,
        cex.lab = 1.25, cex.axis= 1.1, mgp=c(2,0.5,0))
    text(x=100, y=0.05,labels=LOCYR,cex=1.5)
    arrows(S1$DAP, S1$fPAR-S1$fPAR_sd, S1$DAP, S1$fPAR+S1$fPAR_sd, 
           length=0.05, angle=90, code=3)
    points(x = S1$DAP, y = S1$fPAR, pch = 2)
    colLEG <- "black"; colLWD <- 1.5; colLTY <- 1; colPCH <- 2

  if(length(S2)>1){
    lines(S2$DAP, S2$fPAR, col = red, lwd = 1.5, lty = 2)
    arrows(S2$DAP+3, S2$fPAR-S2$fPAR_sd, S2$DAP+3, S2$fPAR+S2$fPAR_sd, 
           length=0.05, angle=90, code=3, col = red)
    points(S1$DAP, S1$fPAR, pch = 5)
    colLEG <- c(colLEG, red); colLWD <- c(colLWD, 1.5); colLTY <- c(colLTY, 2); colPCH <- c(colPCH, 5)
  }
  if(length(S3)>1){
      lines(S3$DAP, S3$fPAR, col=blue, lwd=1.5, lty=4)
      arrows(S3$DAP-3, S3$fPAR-S3$fPAR_sd, S3$DAP-3, S3$fPAR+S3$fPAR_sd, 
             length=0.05, angle=90, code=3,col = blue)
      points(S3$DAP, S3$fPAR, col = blue, pch = 19)
      colLEG <- c(colLEG, blue); colLWD <- c(colLWD, 1.5); colLTY <- c(colLTY, 4); colPCH <- c(colPCH, 19)
  }    
  if(LEG){
    legend("bottomright", legend=LEGTEXT, col = colLEG, lwd = colLWD, 
            lty = colLTY, pch = colPCH, bty = "n")
  }
}

DATA <- read.csv(lcp("ex/fPAR_Benue_Edo.csv"))

par(mfrow=c(2,2), oma=c(3,3,1,1), mar=c(1.5,1.5,0,0))
DAP_vs_IPAR(DATA, "Edo", 2016, TRUE)
DAP_vs_IPAR(DATA, "Benue", 2016)
DAP_vs_IPAR(DATA, "Edo", 2017)
DAP_vs_IPAR(DATA, "Benue", 2017)
mtext("Days after planting", side = 1, line = 1, outer = TRUE)
mtext("IPAR, f", side = 2, line = 1, outer = TRUE, at = 0.25)
mtext("IPAR, f", side = 2, line = 1, outer = TRUE, at = 0.75)
```


### Figure 3

Using data on cumulative intercepted PAR, combining measure wheater data and interpolated measurements of LI from AccurPAR data.


First the regression model


```{r anova_fig3}
Bio_cumIPAR <- read.csv(lcp("ex/Bio_cumIPAR_NFPFKF_BenueEdo.csv"))

i <- which(Bio_cumIPAR$Loc_yr == "Edo_2016")
ModBE <- lm(BiomassDM.gm.2 ~ -1 + Cum_PAR.MJm.2, data=Bio_cumIPAR[i,])
summary(ModBE)
anova(ModBE)
```

And now the plot 

```{r fig3, fig.width=7, fig.height=7}

Fig_3 <- function(bp, model){
  
  par(oma=c(2,2,2,2),mar=c(5,5,0,0))
  bp$Loc_yr <- factor(bp$Loc_yr)
  #par(mfrow=c(1,2))
  i <- which(bp[,"Loc_yr"]=="Edo_2016" )
  plot(bp[i,"Cum_PAR.MJm.2"], bp[i,"BiomassDM.gm.2"], 
       xlab = "Cumulative IPAR (MJ/m²)", ylab = "Biomass (g/m²)",
       cex.lab=1.3, cex = 1.0, pch =1, xlim = c(0,2600), ylim = c(0,9000))
  i <- which(bp[,"Loc_yr"]=="Edo_2017" )
  points(bp[i,"Cum_PAR.MJm.2"], bp[i,"BiomassDM.gm.2"], pch=16)
  
  i <- which(bp[,"Loc_yr"]=="Benue_2016" )
  points(bp[i,"Cum_PAR.MJm.2"], bp[i,"BiomassDM.gm.2"], pch =2)
  i <- which(bp[,"Loc_yr"]=="Benue_2016" )
  points(bp[i,"Cum_PAR.MJm.2"], bp[i,"BiomassDM.gm.2"], pch=17)
  legend('bottomright', c('Edo, 2016', 'Edo, 2017', 'Benue, 2016', 'Benue, 2017'), 
         pch = c(1,16,2,17), bty = 'n')
  
  abline(model)
  legend('topleft', c('Edo 2016', 'R² = 0.98, P< .001', 'y = 2.76x'), bty = 'n')
}


Fig_3(Bio_cumIPAR, ModBE)

```


### Figure 4

```{r fig4, fig.width=10, fig.height=5.5}
#Data output from LINTUL-cASSAVA simulations
#Analysing Simulated and observed yield 
#Read in the data needed
#plant parts g DM/m2
BIOMASS <- read.csv(lcp('ex/Biomass_NPKuptake.csv'))
BIOMASS <- BIOMASS[BIOMASS$Treatment == "NfPfKf", ]
BIOM_Edo16 <- BIOMASS[which(BIOMASS$Location == 'Edo' & BIOMASS$Year == 2016), ]
BIOM_Edo17 <- BIOMASS[which(BIOMASS$Location == 'Edo' & BIOMASS$Year == 2017),]
BIOM_BEN16 <- BIOMASS[which(BIOMASS$Location == 'Benue' & BIOMASS$Year == 2016),]
BIOM_BEN17 <- BIOMASS[which(BIOMASS$Location == 'Benue' & BIOMASS$Year == 2017),]
BIOM_CRS16 <- BIOMASS[which(BIOMASS$Location == 'CRS' & BIOMASS$Year == 2016),]
BIOM_CRS17 <- BIOMASS[which(BIOMASS$Location == 'CRS' & BIOMASS$Year == 2017),]

# leaves

LeavesStemRootDynamics <- function(Series1, Series2=NULL,
                    BIOMASS=NA, TITLE, LEG=FALSE, 
                    LEGTEXT=c('Potential', 'Water lim.', "Observed")){

  red <- rgb(1, 0, 0, alpha = 0.8)

  plot(Series1$time, Series1$WLVG, type="l", xlab = "time (days)",
       ylab = expression('Green leaves, g DM m'^-2), ylim = c(0,500),
       lwd = 2, lty = 1, cex.lab  = 1.25, cex.axis= 1.1, mgp=c(2,0.5,0))
  text(x=median(Series1$time), y=490,labels=TITLE,cex=1.5)
  colLEG <- "black"
  colLWD <- 2
  colLTY <- 1
  colPCH <- NA
  if(length(Series2) > 1){
    lines(Series2$time, Series2$WLVG, col = red, lwd = 2, lty = 2)
    colLEG <- c(colLEG, red)
    colLWD <- c(colLWD, 2)
    colLTY <- c(colLTY, 2)
    colPCH <- c(colPCH, NA)
  }
  if(length(BIOMASS)>1){
    arrows(BIOMASS$time, BIOMASS$DMLeaves.g.DM.m2 - BIOMASS$std_DMLeaves.g.DM.m2, 
           BIOMASS$time, BIOMASS$DMLeaves.g.DM.m2 + BIOMASS$std_DMLeaves.g.DM.m2, 
           length=0.05, angle=90, code=3)
    points(BIOMASS$time, BIOMASS$DMLeaves.g.DM.m2, lwd = 2.5, pch = 19)
    colLEG <- c(colLEG, "black")
    colLWD <- c(colLWD, 2.5)
    colLTY <- c(colLTY, NA)
    colPCH <- c(colPCH, 19)
  }
  
  ##############################################
  
  #STEMS
  
  #par(mfrow=c(1,1), oma = c(4,1,1,1), mar = c(0.6,5,0,0))
  plot(Series1$time, Series1$WST, type="l", xlab = "time (days)", 
        ylab = expression('Stems, g DM m'^-2), ylim = c(0,2000), 
        lwd = 2, cex.lab  = 1.25, cex.axis= 1.1, mgp=c(2,0.5,0))
  
  if(length(Series2)>1){
    lines(Series2$time, Series2$WST, col=red, lwd=2, lty=2)
  }
  if(length(BIOMASS)>1){
    arrows(BIOMASS$time, BIOMASS$DMstems.g.DM.m2-BIOMASS$std_DMstems.g.DM.m2, 
           BIOMASS$time, BIOMASS$DMstems.g.DM.m2+BIOMASS$std_DMstems.g.DM.m2, 
           length=0.05, angle=90, code=3)
    points(BIOMASS$time, BIOMASS$DMstems.g.DM.m2, lwd = 2.5, pch = 19)
  }
  if(LEG){
    legend("topleft", legend=LEGTEXT, col = colLEG, lwd = colLWD, 
           lty = colLTY, pch = colPCH, bty = "n")
  }
  
  ###################################################################
  
  #Storage roots
  
  plot(Series1$time, Series1$WSO, type="l", xlab = "time (days)",
       ylab = expression('Storage roots, g DM m'^-2), ylim = c(0,4000), 
       lwd = 2, cex.lab  = 1.25, cex.axis= 1.1, mgp=c(2,0.5,0))
  if(length(Series2)>1){
    lines(Series2$time, Series2$WSO, col = red, lwd = 2, lty = 2)
  }
  if(length(BIOMASS)>1){
    arrows(BIOMASS$time, BIOMASS$DMRoots.g.DM.m2 - BIOMASS$std_DMRoots.g.DM.m2, 
           BIOMASS$time, BIOMASS$DMRoots.g.DM.m2 + BIOMASS$std_DMRoots.g.DM.m2, 
           length=0.05, angle=90, code=3)
    points(BIOMASS$time, BIOMASS$DMRoots.g.DM.m2, lwd = 2.5, pch = 19)
  }
}


Figure4 <- function(){
  par(mfcol=c(3,5),oma=c(3,3,1,1),mar=c(1.5,1.5,0,0))
  LeavesStemRootDynamics(mod$pot_Edo_2017, mod$wlim_Edo_2017,
                         BIOMASS=BIOM_Edo17, TITLE="Edo '17", LEG=TRUE)
  LeavesStemRootDynamics(mod$pot_CRS_2016, mod$wlim_CRS_2016,
                         BIOMASS=BIOM_CRS16,TITLE="CRS '16")
  LeavesStemRootDynamics(mod$pot_CRS_2017, mod$Wlim_CRS17,
                         BIOMASS=BIOM_CRS17, TITLE="CRS '17")
  LeavesStemRootDynamics(mod$pot_Benue_2016, mod$wlim_Benue_2016,
                         BIOMASS=BIOM_BEN16, TITLE="BEN '16")
  LeavesStemRootDynamics(mod$pot_Benue_2017, mod$wlim_Benue_2017,
                         BIOMASS=BIOM_BEN17, TITLE="BEN '17")
  
  mtext("Days after planting", side = 1, line = 1, outer = TRUE, at = NA)
  mtext("Roots, gDM/m2", side = 2, line = 1, outer = TRUE, at = 0.2)
  mtext("Stems, gDM/m2", side = 2, line = 1, outer = TRUE, at = 0.52)
  mtext("Leaves, gDM/m2 ", side = 2, line = 1, outer = TRUE, at = 0.86)
}

mod <- lapply(mod, \(m) {
    m$time <- as.integer(format(m$date[1], "%j")) + (1:nrow(m)-1)
    m
})


Figure4()
```


### Figure 5

Plotting simulated versus observed storage root yield NfPfKf for all years in Benue, Cross River and Edo  (Fig. 5 in thesis chapter 3)

```{r fig5, fig.width=7, fig.height=7}
#Get the data
measured <- read.csv(lcp("ex/Biomass_NPKuptake.csv"))

measured <- measured[measured$Treatment == "NfPfKf", c("Year","Location","time","DMRoots.g.DM.m2")]
measured[,"DMRoots.t.ha"] <- measured[,"DMRoots.g.DM.m2"] * 0.01 #1E4*1E-6

mod <- lapply(mod, \(m) {
    m$WSOTHA <- m$WSO / 100
    m$Year <- as.integer(format(m$date[1], "%Y"))
    m
})

modelled <- do.call(rbind, mod[grep("^wlim", names(mod))])[, c("Location", "Year", "time", "WSOTHA")]


dc <- merge(measured, modelled, by=c("time", "Year", "Location"))

fit <- lm(WSOTHA ~ -1 + DMRoots.t.ha, data=dc)
sfit = summary(fit)
RMSE=sqrt(mean(sfit$residuals^2))
anova(fit)
print(paste0("Fitted linear regression. Slope =", sfit$coefficients[1], ", 
             R2= ", sfit$r.squared ,", RMSE =", RMSE))


Figure_5 <- function(d){
  i <- which(d$Year==2017 & d$Location=="Benue")
  plot(x=d[i,"DMRoots.t.ha"],y=d[i,"WSOTHA"],
       xlim = c(0,40), ylim = c(0,40), cex.lab=1.5, cex = 1.5,
       xlab = 'Observed storage roots (t DM/ha)',
       ylab = 'Simulated storage roots (t DM/ha)', pch=15)
  i <- which(d$Year==2016 & d$Location=="Benue")
  points(d[i,"DMRoots.t.ha"], d[i, "WSOTHA"], pch=0, cex=1.5)
  i <- which(d$Year==2017 & d$Location=="Edo")
  points(d[i,"DMRoots.t.ha"], d[i, "WSOTHA"], pch=16,cex=1.5)
  i <- which(d$Year==2016 & d$Location=="CRS")
  points(d[i,"DMRoots.t.ha"], d[i, "WSOTHA"], pch=2, cex=1.5)
  i <- which(d$Year==2017 & d$Location=="CRS")
  points(d[i,"DMRoots.t.ha"], d[i, "WSOTHA"], pch=17, cex=1.5)
  lines(c(-100,100), c(-100,100), lwd=2)
  lines(c(-100,100), c(-100,100)*sfit$coefficients[1], col="red",lty=3,lwd=3)
  
  legend("bottomright", c("Edo, 2017", "CRS, 2016", "CRS, 2017", "Benue, 2016" ,"Benue, 2017"),
         cex=1.2, pch = c(16, 2,17,0,15), bty="n") 
  #text(x=10,y=35,labels=paste0("slope = ",round(sfit$coefficients[1],2),", RMSE=",round(RMSE,2),", R²=",round(sfit$r.squared,2)))
  legend('topleft', c( paste0('y = ', round(sfit$coefficients[1],2),'x'),
                       paste0('R² = ', round(sfit$r.squared,2)),
                       paste0('RMSE = ',round(RMSE,2))), cex=1.2, bty='n')
}

Figure_5(dc)

```