04-soilwater.Rmd

# Soil Water 

## Aim of this document

1. Explore the soil water data  
2. Try if multiple change point analysis could help to estimate DULs or CLLs
3. Try if `nls` can help estimate the KLs

## Experiment details 

The experiment detailed in [@Sim2014] thesis. Here is a short description to 
note.

a. **2** years experiment.
b. **2** sites
   a. Lincoln Uni
   b. Ashely Dene
c. **10** sowing date in both sites \@ref(tab:sowingdates). 
   a. five sowing date in the first year
   b. five sowing date in second year
   
|Year 1 (Seedling crop) | Year 2 (Seedling & Regrowth)|
|--------| ------------|
|SD1, SD2, SD3, SD4, SD5| The first five sowing date Regrowth|
||SD6, SD7, SD8, SD9, SD10|The second five sowing date seedling|
_Note:_ The setup helps to compare the seedling and regrowth in the same conditions in Year 2. 

d. **1** cultivar.
e. **4** reps each site. 
f. **2** soil types
   a. stone free - plant available water capaciy (PAWC) 360mm/2.3m
   b. ~~stony - PAWC 240/2.3m~~ This type is for a grazing trial. probably not a good idea for a starting?
   c. very stony 130mm/2.3m 
g. **22** layers of soil water measurements - 10 cm interval below top 20 cm
   a. top 20 cm was measured by TDR. 
   c. remaining 21 layers was meastured by Neutron probe
h. **1** follow treat in both sites - absolute soil evaporation 


```{r import pkg, include=FALSE}
# knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE,
#                       fig.align = "center", fig.width = 10, fig.height = 7, 
#                       dpi = 300)
options(scipen = 999)
#Constants
outlier.colour = "#ff0000"
source(here::here("02Scripts/R/packages.R"))
source(here::here("02Scripts/R/functions.R"))


```

```{r read sim, warning=FALSE}
id_vars <- c("Experiment", "SowingDate", "Clock.Today", "DAS")
con <- dbConnect(SQLite(), here::here("01Data/ProcessedData/Richard.sqlite3"))

data_sw <- dbReadTable(con, "SoilWater") %>% 
  as.data.table()
data_sw[, Clock.Today := as.Date(Clock.Today,  origin = "1970-01-01")]
SowingDates <- dbReadTable(con, "SowingDates") %>% 
  as.data.table()
SowingDates[, Clock.Today := as.Date(Clock.Today,  origin = "1970-01-01")]

value_vars <-  grep("SWmm\\.\\d.", colnames(data_sw), value = TRUE)
sw_mean_new <- colwise_meanSW(data_sw)

```
```{r sowingdates}
knitr::kable(
  list(SowingDates[1:10],
       SowingDates[11:20]),booktabs = TRUE,
  caption = 'Sowing dates details in two sites.') %>% 
    kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
```

## EDA 

First Layer

```{r scatter, echo=FALSE, message=FALSE, warning=FALSE}
sw_mean_new %>% 
  ggplot(aes(Clock.Today, SWmm.1..mean)) +
  geom_point()+ 
  facet_wrap(SowingDate ~ Experiment)
```


```{r hist, echo=FALSE, message=FALSE, warning=FALSE}
sw_mean_new %>% 
  ggplot(aes(SWmm.1..mean)) + 
  geom_histogram(binwidth = 3) + 
  facet_wrap( ~ Experiment)
```

Figure \@ref(fig:relativeSW) shows the relative SW measurements in each layer and scatter plot. 

```{r relativeSW,fig.cap="Soil water measurements in each layer in relation to the layer DULs",echo=FALSE, message=FALSE, warning=FALSE}
VWC <- grep("VWC", x = colnames(sw_mean_new), value = TRUE)
VWCcols <- c(id_vars, VWC)
DT_VWC <- sw_mean_new[,..VWCcols] %>% 
  melt.data.table(id.vars = id_vars, 
                  value.name = "SW",
                  variable.name = "Depth",
                  variable.factor = FALSE)
DT_VWC[, Depth:= as.integer(gsub("\\D","",Depth))]
DT_VWC[, DUL := max(SW), by = .(Experiment, SowingDate, Depth)
       ][, relativeSW:= SW/DUL]

DT_VWC %>% 
  ggplot(aes(Clock.Today, relativeSW, color = SowingDate)) +
  geom_point() +
  geom_line()+
  scale_x_date(date_breaks = "8 weeks") +
  geom_smooth()

```

**Subset** the autumn rain.  

```{r ADonly, fig.cap="Ashley Dene, soil water measurements in each layer in relation to the layer DULs",echo=FALSE, message=FALSE, warning=FALSE}
AD <- DT_VWC[ Experiment == "AshleyDene" &
          Clock.Today %between% c("2011-04-01", "2012-04-30")]

# AD[order(SW, decreasing = TRUE), .SD[1] , by = c(id_vars, "Depth")]


AD %>% 
  ggplot(aes(Clock.Today, relativeSW, group = Depth)) +
  geom_point() +
  scale_x_date(date_breaks = "8 weeks") +
  geom_smooth(formula = 'y~x', span = 0.25, se = FALSE, method = 'loess')
```


```{r ADonlydensity, fig.cap="Ashley Dene, hist relative SW",echo=FALSE, message=FALSE, warning=FALSE}
AD %>% 
  ggplot(aes(relativeSW)) +
  geom_histogram(bins = 50)+
  geom_density()
```

## Multiple change point analysis 

### Try a simple model 

**Define a model for multiple change points analysis**

**Red** line is the mean.  
**Black** line is the measured value at the closet date to obtain the mean value. 

```{r try mcp, message=FALSE, warning=FALSE, include=FALSE}
model = list(
  relativeSW ~ 1 + sigma(1),  # plateau (int_1)
  ~ 0 + DAS ,       # joined slope (time_2) at cp_1, could be a exp decay function.
  ~ 0

)
# Define the prior
prior = list(
  # Intercept should be between 0.5 and 1; less than 1
  ## Evidence in the hist above
  int_1 = "dnorm(0.7, 1) T(, 1)" 
)
dt <- AD[SowingDate == "SD1"]
# Fit it. 
fit = mcp(model, data = dt, cores = 3, prior = prior)
    # Extract the cp one 
cp_1_est = as.data.table(mcp::fixef(fit))
setkey(dt, DAS)
close_das = dt[dt[J(cp_1_est$mean[1]), roll = 'nearest', which = TRUE]]$DAS

P = plot(fit) +
      geom_vline(xintercept = c(cp_1_est$mean[1],close_das), 
                 color =  c("red","black")) +
      # ggtitle(paste0("Layer ",1))+
  theme_water()
```


```{r mcpplotrs, echo=FALSE, message=FALSE, warning=FALSE}
P
```



```{r mcprs, echo=FALSE, message=FALSE, warning=FALSE}
cp_1_est %>% 
  kable(caption = "Asheley Dene SD1 first attempt results.") %>% 
    kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
```
```{r tryALayer2 , echo=FALSE, message=FALSE, warning=FALSE, include=FALSE}
dt <- AD[SowingDate == "SD1" & Depth == 2]
# Fit it. 
fit = mcp(model, data = dt, cores = 3, prior = prior)
    # Extract the cp one 
cp_1_est = as.data.table(fixef(fit))
setkey(dt, DAS)
close_das = dt[dt[J(cp_1_est$mean[1]), roll = 'nearest', which = TRUE]]$DAS

P = plot(fit) +
      geom_vline(xintercept = c(cp_1_est$mean[1],close_das), 
                 color =  c("red","black")) +
      ggtitle(paste0("Layer ",2))+
  theme_water()
```


```{r layer2plot , echo=FALSE, message=FALSE, warning=FALSE}
P
```

```{r trylayer2rs, echo=FALSE, message=FALSE, warning=FALSE}
cp_1_est %>% 
  kable(caption = "Asheley Dene SD1 layer 2 first attempt results.") %>% 
    kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
```

```{r verfiylayer13AD , echo=FALSE, message=FALSE, warning=FALSE, include=FALSE}
dt <- AD[SowingDate == "SD1" & Depth == 13]
# Fit it. 
fit = mcp(model, data = dt, cores = 3, prior = prior)
    # Extract the cp one 
cp_1_est = as.data.table(fixef(fit))
setkey(dt, DAS)
close_das = dt[dt[J(cp_1_est$mean[1]), roll = 'nearest', which = TRUE]]$DAS

P = plot(fit) +
      geom_vline(xintercept = c(cp_1_est$mean[1],close_das), 
                 color =  c("red","black")) +
      ggtitle(paste0("Layer ",3))+
  theme_water()
```


```{r verfiylayer13ADplot , echo=FALSE, message=FALSE, warning=FALSE, include=FALSE}
P
```


```{r verfiylayer13ADrs, echo=FALSE, message=FALSE, warning=FALSE}
cp_1_est %>% 
  kable(caption = "Asheley Dene SD1 layer 13 first attempt results.")
```


### Forming a function

```{r}
estimate_DUL <- function(DT, sowingdate, layer, mcpmodel = model){
  dt <- AD[SowingDate == sowingdate & Depth == layer]
  # Fit it. 
  fit = mcp(mcpmodel, data = dt, cores = 3, prior = prior)
  # Extract the cp one 
  cp_1_est = as.data.table(fixef(fit))
  setkey(dt, DAS)
  close_das = dt[dt[J(cp_1_est$mean[1]), roll = 'nearest', which = TRUE]
                 ]
  l <- list(fit,cp_1_est, close_das)
  names(l) <- c("model","cp1", "nearDAS")
  return(l)
}

```

### Test the function

```{r testfun, echo=FALSE, message=FALSE, warning=FALSE, include=FALSE}
ADSD1 <- estimate_DUL(AD, "SD1", 1, mcpmodel = model)
cp1_mean <- ADSD1$cp1$mean[1]
cp1_nearDAS <- ADSD1$nearDAS$DAS
cp1_int1 <- ADSD1$cp1$mean[4]
```


```{r testfunlayer1, fig.cap="Visualise the mcp results.", echo=FALSE, message=FALSE, warning=FALSE}
plot(ADSD1$model) +
  geom_vline(xintercept = c(cp1_mean,cp1_nearDAS), 
             color =  c("red","black")) +
  geom_hline(yintercept = cp1_int1, color = "blue") +
  annotate(geom = "text", 
           x = cp1_mean * 1.03, 
           y = cp1_int1 * 1.02, label = round(cp1_int1, digits = 3)) +
  ggtitle(paste0("Layer ",1))+
  theme_water()
```

### Make a loop


Go through all layers and estimate a DUL.  

```{r, include=FALSE, echo=FALSE, message=FALSE, warning=FALSE}
list_AD <- AD[, list(data=list(.SD)), by = .(Experiment, SowingDate, Depth)]
# list_AD$data[1]

list_AD_models <- list_AD[, esti := list(apply(.SD, 1, function(x){
  l <- estimate_DUL(x[["data"]], x[["SowingDate"]], x[["Depth"]])
  return(l)
})), by = .(Experiment)]
# list_AD[, .SD[1]]

plot(list_AD_models$esti[[1]]$model)

TEST <- list_AD_models[,.SD[1]][, graphs := list(lapply(esti, function(x){
  cp1_mean <- x$cp1$mean[1]
  cp1_nearDAS <- x$nearDAS$DAS
  cp1_int1 <- x$cp1$mean[4]
  depth <- x$nearDAS$Depth
  p <- plot(x$model) +
    geom_vline(xintercept = c(cp1_mean,cp1_nearDAS), 
               color =  c("red","black")) +
    geom_hline(yintercept = cp1_int1, color = "blue") +
    annotate(geom = "text", 
             x = cp1_mean * 1.03, 
             y = cp1_int1 * 1.02, label = round(cp1_int1, digits = 3)) +
    ggtitle(paste0("Layer ", depth))+
    theme_water()
  return(p)
})) ]

```


```{r, echo=FALSE, include=FALSE}
list_AD_models_graphs <- list_AD_models[, graphs := list(lapply(esti, function(x){
 
  cp1_mean <- x$cp1$mean[1]
  cp1_nearDAS <- x$nearDAS$DAS
  cp1_int1 <- x$cp1$mean[4]
  depth <- x$nearDAS$Depth
  sd <- x$nearDAS$SowingDate

   # cat("Process layer", depth, "\r\n")
  p <- plot(x$model) +
    geom_vline(xintercept = c(cp1_mean,cp1_nearDAS), 
               color =  c("red","black")) +
    geom_hline(yintercept = cp1_int1, color = "blue") +
    annotate(geom = "text", 
             x = cp1_mean * 1.03, 
             y = cp1_int1 * 1.02, label = round(cp1_int1, digits = 3)) +
    ggtitle(paste0("Layer ", depth, "-", sd))+
    theme_water()
  return(p)
})) ]
```

```{r, eval=FALSE}
list_AD_models_graphs$esti[[5]]$nearDAS
list_AD_models_graphs$esti[[5]]$cp1
list_AD_models_graphs$graphs[[5]]
```

## What are the estimated DULs?

```{r posteriroDULs, echo=FALSE, include=FALSE}
list_AD_models[, results:= lapply(esti, function(x){
  cp1_int1 <- x$cp1$mean[4]
  DT = x$nearDAS[, esti_DUL := cp1_int1]
  return(DT)
  })]

estimates <- rbindlist(list_AD_models$results)
```
```{r estimatedDULs}
estimates[, posterior_DUL := DUL * esti_DUL] %>% 
  kable(caption = "Estimated DUL in VWC.") %>% 
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")

```

### How much different it made? 

```{r comparisonduls, fig.cap="A comparison of three methods: average the max, mcp estimation and arbitrary 0.95 for estimating DULs", echo=FALSE, warning=FALSE, message=FALSE}
DUL_shift <- estimates %>% 
  ggplot(aes(Depth)) +
  geom_point(aes(y = DUL, color = "DUL")) + 
  geom_smooth(aes(y = DUL, color = "DUL"), se = FALSE, method = 'loess', formula = 'y ~ x') +
  geom_point(aes(y = posterior_DUL, color = "POSTDUL")) + 
  geom_smooth(aes(y = posterior_DUL, color = "POSTDUL"), se = FALSE, method = 'loess', formula = 'y ~ x') +
   geom_point(aes(y = DUL * 0.95, color = "arbitraryDUL")) + 
  geom_smooth(aes(y = DUL * 0.95, color = "arbitraryDUL"), se = FALSE, method = 'loess', formula = 'y ~ x') +
  coord_flip()+
  scale_x_reverse(limits = c(23, 0), expand = c(0,0)) 
  
DUL_shift +
  theme_water()+
  scale_color_manual(values = c("DUL" = "#FF0000", "POSTDUL" = "#0000FF", 
                                "arbitraryDUL" = "#000100"))
```

### Output all graphs for inspections

```{r, eval=FALSE, echo=TRUE}
lapply(list_AD_models$graphs, function(x){
   name <- unique(x$data[,c("Experiment", "SowingDate", "Depth")]) %>%
    unlist() %>%
    paste(., collapse = "-")
  ggsave(filename = here::here("<path>", paste0(name, ".png")), plot = x, 
         dpi = 300)
})
```


### Model2 

May or not worth to try a second model? 


```{r model2, echo=TRUE}
model_2 = list(
  relativeSW ~ 0 + DAS,  # plateau (int_1)
  ~ 1  ,       # joined slope (time_2) at cp_1, could be a exp decay function.
  ~ 0 + DAS
  
)
prior_2 = list(
  # Intercept should be between 0.5 and 1; less than 1
  ## Evidence in the hist above
  int_2 = "dnorm(0.7, 1) T(, 1)" 
)
```


## The final version workflow for DUL, LL and SW via average max and min. 

**DUL AND LL CAN BE SOWINGDATES DEPEND - apsim is a point model.  **

```{r outuputDULLL}
id_vars <- colnames(data_sw)[!colnames(data_sw) %in% grep("SWmm\\.\\d.+", 
                                                          colnames(data_sw), 
                                                          value = TRUE)]
value_vars <- colnames(data_sw)[!colnames(data_sw) %in% id_vars]
id_vars <-  c("Experiment","Clock.Today", "Season","SowingDate", "DAS")
DUL_LL <- data_sw[, lapply(.SD, mean, na.rm = TRUE), 
                  by = id_vars, .SDcols = value_vars # average 4 reps
                  ][, unlist(lapply(.SD, max_min), recursive = FALSE), 
                    by = .(Experiment, SowingDate), .SDcols = value_vars]

melted_DUL_LL <- melt(DUL_LL, 
                      id.var = c("Experiment","SowingDate"), 
                      variable.factor = FALSE) 
DUL_LL_SDs <- melted_DUL_LL[, (c("Depth", "variable")) := tstrsplit(variable, "\\.{2}")] %>% 
  dcast(Experiment +  SowingDate + Depth ~ variable)


# Sort out the sowing date 

SW_mean <- data_sw[, lapply(.SD, mean, na.rm = TRUE), by = id_vars, .SDcols = value_vars]
SW_initials <- SW_mean[SowingDates, on = c("Experiment", "SowingDate", "Clock.Today"), roll = "nearest"]
SW_initials <- melt(SW_initials, 
                    id.vars = id_vars, variable.factor = FALSE,
                    variable.name = "Depth", value.name = "PAWC"
                    )
# Clean the depth column
SW_initials[, Depth := as.integer(gsub("\\D", "", Depth))]
DUL_LL_SDs[, Depth := as.integer(gsub("\\D", "", Depth))]
SW_DUL_LL <- merge.data.table(SW_initials, DUL_LL_SDs, by = c("Experiment", "SowingDate", "Depth"),
                              all.x = TRUE)

# Check if anything unusual 
# SW_DUL_LL[DUL<PAWC]
# SW_DUL_LL[LL>PAWC]

# Order the layers
SW_DUL_LL[, ':='(Depth = as.integer(gsub("\\D", "", Depth)),
                 SowingDate = forcats::fct_relevel(as.factor(SowingDate), paste0("SD", 1:10)))]
# order the order - CRITICAL 
setkey(SW_DUL_LL, Experiment, SowingDate, Depth)
kable(SW_DUL_LL, caption = "Initial soil parameters for 22 layers in 10 different sowing dates.")%>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "500px")
#Make a version that can be plugged into apsimx 
DUL_LL_SDsVWC <- copy(DUL_LL_SDs)[, ':='(DUL = round(DUL/100, digits = 3),
                                         LL = round(LL/100, digits = 3))
                            ]
```

try to modify apsimx file in place. **Edit** feature. 
_Note_: apsimx is not working properly in the bash env.  


```{r}
data_sw[,.(Experiment, Clock.Today, SWmm.1.)] %>% 
  ggplot(aes(Clock.Today, SWmm.1.)) +
  geom_point() +
  facet_grid(~ Experiment)
```

### Examine the DUL and LL with original observations 



```{r manualcheckboundaries, echo=FALSE, fig.cap="Manual visualisation of boundaries of PAWC for each sowing date.", warning=FALSE, message=FALSE}
disered_cols = c(id_vars, value_vars)
water_long = melt(data_sw[,..disered_cols],
                  id.var = id_vars, variable.name = "Depth", 
                  variable.factor = FALSE) 
water_long[, Depth := as.integer(gsub("\\D","", Depth))
           ]
long = merge(water_long, SW_DUL_LL,
      by = c("Experiment", "Clock.Today", "SowingDate", "Season", "DAS", "Depth"), all= TRUE)

long %>% 
  ggplot(aes(Clock.Today, value))+ 
  geom_point(color = "grey") +
  geom_hline(aes(yintercept = DUL), colour = "blue", alpha = 0.8) +
  geom_hline(aes(yintercept = LL), colour = "red", alpha = 0.8) +
  facet_wrap(~Experiment + SowingDate) +
  theme_water()
```
looks reasonable. but is this too busy. 


```{r histDULLL, echo=FALSE, fig.cap="Hist for DUL and CLL. Overlapping suggest some crops yet to complete drying cycle."}
SW_DUL_LL %>% 
  ggplot()+ 
  geom_histogram(aes(DUL, fill = "DUL"),bins = 30)+
  geom_histogram(aes(LL, fill = "LL"), alpha = 0.5,bins = 30) +
  facet_wrap(~Experiment ) +
  theme_water()

```

```{r swinitial, fig.cap="Initial soil water content.", echo=FALSE}
intial_long = melt(SW_DUL_LL, id.vars = c(id_vars, "Depth"), variable.factor = FALSE)

intial_long %>% 
  ggplot()+ 
  geom_point(aes(Depth, value, color = SowingDate), size = 3) +
  geom_line(aes(Depth, value,group = SowingDate, color = SowingDate)) +
  facet_wrap(~Experiment + variable) +
  coord_flip()+
  scale_x_reverse() +
  theme_water()
```


## Estimate KLs

Started at 20200531

Visulise the data to have a feel.

```{r , fig.height=7, fig.width=10, fig.cap="Selected SW in selected SD in Selected site.", echo=FALSE}
SW_mean <- data_sw[, lapply(.SD, mean, na.rm = TRUE), by = id_vars, .SDcols = value_vars]

SD <- "SD2"
Site <- "AshleyDene"
SW_mean[SowingDate == SD & Experiment == Site] %>% 
  ggplot(aes(Clock.Today, SWmm.1.)) +
  geom_point() +
  geom_smooth(span = 0.25,method = 'loess' , formula = 'y ~ x') +
  scale_x_date(breaks = "2 weeks") +
  theme_classic() +
  geom_text(aes(label = Clock.Today), check_overlap = TRUE) +
  theme(axis.text.x = element_text(angle = 90))



```

**How to deal with outlines?!!**

Do I just assume that all the water was sucked up by the
`Areas within 1m of the neutron probe access tubes were also avoided.`
no harvest around NP? always full canopy?

https://douglas-watson.github.io/post/2018-09_exponential_curve_fitting/

Eliminate the noise data   
1. full canopy cover - LAI > 1  
2. No rainfall during the period. 

```{r,warning=FALSE}
dt1 <- SW_mean[SowingDate == "SD1" & Experiment == "AshleyDene" & Clock.Today > "2011-07-27",
               .(Clock.Today, SW = SWmm.1.)]
qplot(x = Clock.Today, y = SW, data = dt1)
# 2011-07-27 Had no rainfall, the SW was the same like the previous measurement. 
# Therefore, a good starting point. 

```
sw
need to ask Hamish about the model-fitting exercises. [@Brown2009]


```{r, fig.width=10, fig.height=7, echo=FALSE, fig.cap="Label the dates for fitting the curves.",warning=FALSE}
con <- dbConnect(SQLite(), here::here("01Data/ProcessedData/Richard.sqlite3"))
biomass <- dbReadTable(con, "biomass") %>% 
  as.data.table()
met_AD <- dbReadTable(con, "met_AshleyDene") %>% 
  as.data.table()
dbDisconnect(con)
SD1fullCano <- biomass[,.(Experiment, SowingDate, 
                          Clock.Today = as.Date(Clock.Today), 
                          Rotation = as.character(Rotation.No.),Harvest.No., LAI..FW.)
                       ][Experiment == "AshleyDene"
                         ][LAI..FW.>1 & LAI..FW. < 1.5
                           ][, list(LAI_mean = mean(LAI..FW.)), 
                             by = .(SowingDate, Clock.Today,Rotation,Harvest.No.)]
SD1fullCano %>% 
  ggplot(aes(Clock.Today, LAI_mean, color = Rotation, shape = SowingDate))+
  geom_point()
# The dates that has full canopy clousre for SD1 
dates_interested <- unique(SD1fullCano$Clock.Today)
water_SD1 <- data_sw[Experiment == "AshleyDene"]
data_sw_mean <- water_SD1[, lapply(.SD, mean, na.rm = TRUE), by = id_vars, .SDcols = value_vars]
water_mean <- melt(data_sw_mean, 
                   id.vars = id_vars,
                   variable.factor = FALSE)
# bring in the rainfall 
met_AD <- met_AD[, day := seq(from = 152, length.out = .N, by = 1)
                 ][, Clock.Today := as.Date(day, origin = "2010-01-01")
                   ][, rain := ifelse(rain == 0.0, NA, rain)
                     ][,.(Clock.Today, rain)]

pallete <- c("#0000ff", "#ff0000", "black")
names(pallete) <- c("SW", "LAI_mean", "rain")
datesForCurveFitting <- as.Date(c("2011-10-21", "2012-01-16"))
p <- water_mean[SowingDate == "SD1" & variable %in% paste0("SWmm.", 1:3, ".")] %>% 
  ggplot(aes(x = Clock.Today)) +
  geom_vline(xintercept = datesForCurveFitting) +
  geom_point(aes(y = value, color = "SW"))+
  geom_point(data = SD1fullCano, aes(y = LAI_mean*10, color = "LAI_mean")) +
  geom_point(data = met_AD[Clock.Today %in% dates_interested],
             aes(y = rain, color = "rain")) +
  scale_color_manual(values = pallete) + 
  facet_grid(variable ~ ., scales = "free_y")
p
```


```{r, fig.width=10, fig.height=7, echo=FALSE, fig.cap="Zoom in a bit more.",warning=FALSE}
p + coord_cartesian(xlim = as.Date(c("2011-10-21", "2012-06-16")))
```

Richard conversion of met file had one day off   
152th of 2010 should be 2nd of June  2010 

In the thesis,   
arbitrary dates were selected for extraction initiation and finishing.   
p147: _Extraction began in the 0.25 m soil layer on the 25 October 2011 at both sites._  
p150: _The EFV reached the maximum depth of measurement (2.25 m) on 23 December 2011 at Ashley Dene and 11 March 2012 at Iversen 12_
_Crops at both Ashley Dene and Iversen reached maximum extraction depth in late February 2012._
Assumption based on the graph below. soil water content was averaged across
different sowing dates in the second season. Water extraction was observed to
start in early August 2011 for regrowth crops, however a large rainfall event
returned SWC to DUL on **21 October 2011**. Models were fitted from this point
to maximum SWD which occurred on **16 January 2012 and 25 May 2012** at Ashley
Dene and Iversen 12, respectively.

p136   $$SWC_t = LL + PAWC exp(-kl(t-tc)Sc)$$

```{r}
artificial <- water_mean[variable %in% paste0("SWmm.", 1,".") 
                         & SowingDate %in% paste0("SD", 1:5) 
                         & Clock.Today >= "2011-10-21" 
                         & Clock.Today <= "2012-03-01"
                         ][, list(mean_SW = mean(value)),
                           by = .(Experiment,  Clock.Today)] 
```

```{r sliceouttheperiod, fig.cap="Measured mean SWC during the period."}
artificial %>% 
  ggplot(aes(x = Clock.Today)) +
  geom_point(aes(y = mean_SW)) +
  geom_text(aes(y = mean_SW, label = mean_SW))

```

```{r}
artificial[, start := as.Date("2010-10-21")
           ][, t := as.numeric(difftime(Clock.Today,start, units = "day"))][]
pawc <- 27
ll <- 6

fit <- nls(mean_SW ~ SSasymp(t, Asym = pawc, R0 = ll, x), data = artificial)

# fit
```


```{r quicklook, fig.cap="The quick plot to examine the fitting results"}
qplot(t, mean_SW, data = broom::augment(fit)) + geom_line(aes(y = .fitted))
```


```{r}
broom::tidy(fit) %>% 
  dplyr::select(term, estimate) %>% 
  tidyr::spread(term, estimate) %>% 
  dplyr::mutate(alpha = exp(x)) %>% 
  kable(caption = "The fitting number") %>% 
  kable_styling()


```

is this manual curve fitting procedure robust? the rainfall definetly affect the SW in the top soils?  
how can we ensure that the decay is due to the plant transipiration rather than evaporation?? 

### More fitting trys 

20200515  
https://rpubs.com/mengxu/exponential-model   

 $$SWC_t = LL + PAWC * exp(-kl(x)Sc)$$
 why not $$SWC_t = DUL - PAWC * exp(-kl(t-tc)Sc)$$
 because the plant take water up rather than storing water. 
 The SW will need to get smaller over time. 
 
 
```{r layer9fitting, fig.cap="Examine the mean observations for layer 9 over first five sowing dates in Ashley Dene."}
artificial <- copy(water_mean)[variable %in% paste0("SWmm.", 9,".") 
                               & SowingDate %in% paste0("SD", 1:5) 
                               & Clock.Today >= "2011-10-21" 
                               & Clock.Today <= "2012-05-15"
                               ][, list(mean_SW = mean(value)/100),
                                 by = .(Experiment,  Clock.Today, SowingDate)] 
artificial %>% 
  ggplot(aes(x = Clock.Today)) +
  geom_point(aes(y = mean_SW)) +
  geom_text(aes(y = mean_SW, label = mean_SW)) +
  facet_wrap(~ Experiment)
```


```{r}
artificial[, start := as.Date("2011-10-21")
           ][, t := as.numeric(difftime(Clock.Today, start, units = "day"))]
# Nested data into a list column
nested <- artificial[, list(model = list(.SD)), 
                     by = .(Experiment, SowingDate), .SDcols = c("Clock.Today","mean_SW", "start", "t")]
# paras_SW <- DUL_LL_SDs[Experiment == "AshleyDene" 
#                        & SowingDate %in% paste0("SD", 1:5)
#                        & Depth %in% paste0("SW(", 1,")") ]
# paras_SW_SD1to5 <- nested[paras_SW, on = c("Experiment", "SowingDate")]

data.df <- nested$model[[1]]
# Select an approximate $\theta$, since theta must be lower than min(y), and greater than zero
ll.0 <- min(data.df$mean_SW) * 0.5  

# Estimate the rest parameters using a linear model
model.0 <- lm(log(mean_SW - ll.0) ~ t, data=data.df)  
pawc.0 <- exp(coef(model.0)[1])
kl.0 <- coef(model.0)[2]

# Starting parameters
start <- list( pawc = pawc.0, kl = kl.0, ll = ll.0)

model <- nls(mean_SW ~ pawc * exp(kl * t) + ll , data = data.df, start = start)

# Plot fitted curve
# plot(data.df$t, data.df$mean_SW)
# lines(data.df$t, predict(model, list(x = data.df$t)), col = 'skyblue', lwd = 3)
```


```{r fitlayer9, fig.cap="Examine fitting for layer 9 over five sowing dates in AshleyDene"}
qplot(t, mean_SW, data = broom::augment(model)) + 
  geom_line(aes(y = .fitted))
```


```{r}
broom::tidy(model) %>% 
  dplyr::select(term, estimate) %>% 
  tidyr::spread(term, estimate) %>% 
  kable(caption = "Layer 9 fitting results.") %>% 
  kable_styling()
 
```

Seems too much effort to figure out the starting point by oneself.? 
`selfstarting` might be the way to go

Asked Rodelyn about the `nls` and `start` point preparation.   
the online document takes two steps for fitting the curve:  
1. prepare the `nls` a starting point by transforming it to a linear function.  
2. apply the origin function in `nls` with the `start` 

```{r, echo=FALSE, include=FALSE}
water_mean <- data_sw[SowingDate %in% paste0("SD", 1:10) 
                    ][, lapply(.SD, mean, na.rm = TRUE), by = id_vars, .SDcols = value_vars]
water_mean <- melt(water_mean, 
                   id.vars = id_vars,
                   variable.factor = FALSE)[, value := value/100]
# Nested data into a list column
nested <- water_mean[Clock.Today >= "2011-10-21" 
                     & Clock.Today <= "2012-05-15"
                     ][, list(data = list(.SD)), 
                       by = .(Experiment, SowingDate, variable), 
                       .SDcols = c("Clock.Today","value")]
# Relevel the order 
nested[, ':='(SowingDate=forcats::fct_relevel(as.factor(SowingDate), paste0("SD", 1:10)),
              variable = forcats::fct_relevel(as.factor(variable), paste0("SWmm", 1:22,".")))]

# Add a t value to represent the days difference
nested <- nested[, data := lapply(data, function(x){
  x <- x[, c(unique(.SD))
         ][, start := as.Date("2011-10-21")
           ][, t := as.numeric(difftime(Clock.Today, start, units = "day"))]
  x
  }
  )]
# paras_list <- vector("list", nrow(nested))
nested[, paras_list:= lapply(data, function(x) {
  dt = x
# for(i in seq_len(nrow(nested))){
  # dt = nested$data[[i]]
  # name = paste(nested$Experiment[[i]], 
  #               nested$SowingDate[[i]],
  #               nested$variable[[i]],sep = "_")
  # cat(name, "\r\n")
  # Lower limit starting point 
  ll.0 =  min(dt$value)/2
  # Use a simple liner model to get the starting points for the decay rate _kl_
  # The liner model has to be log transferred first 
  model.0 = lm(log(value - ll.0) ~ t, data = dt)  
  pawc.0 = exp(coef(model.0)[1]) # The starting point for interception
  kl.0 = coef(model.0)[2] # kl starting point
  #  Starting parameters
  start = list(PAWC = pawc.0, kl = kl.0, LL = ll.0)
  model = tryCatch({
    nls(value ~ PAWC * exp(kl * t) + LL , data = dt, start = start)
    },
    error = function(cond){
      print("Warning!!! nls can't describe part of the data.")
      return(NA)
    })
  # paras_list[[i]] <- as.data.table(broom::tidy(model))
  return(model)

# } # the end of for loop 
  })]
```

```{r tidy model fitting, warning=FALSE, message=FALSE, include=FALSE}
nested$paras_list[[1]] 
nested[, augment := lapply(paras_list, function(x){
  x <- tryCatch(broom::augment(x), 
                error = function(cond){
                  message(cond)
                  return(NA)
                })
  x
  })]
nested[, tidied := lapply(paras_list, broom::tidy)]
# debug a dataset that failed nls 
# "./R/20200515curveFittingKLs.R"
```


```{r fittingrs,out.width='200%', fig.width=14, fig.asp=.75, fig.cap="Fitting results for all layers and all sowing dates"}
kl_present <- nested[!is.na(paras_list)]
unnested_plot <- kl_present[, unlist(augment, recursive = FALSE), 
                            by = .(Experiment, SowingDate, variable)]

p <- unnested_plot %>%
  ggplot(aes(t, value)) +
  geom_point() + 
  geom_line(aes(y = .fitted), color = "#ff0000")+ 
  facet_grid(variable ~ Experiment + SowingDate) +
  ggthemes::theme_few()
p
# qplot(t, mean_SW, data = broom::augment(model)) + geom_line(aes(y = .fitted))
# ggsave("../05figures/SoilWaterEDA/kl_curveFitting.pdf", dpi = 300, height = 10,width = 20)

```

### Curve fitting is not going to work 

kl should be link to the root turnover?? **YES**  
But it can't be dynamic to reflect the real world stuiation in NZ since the cold temperature could limit the crop water uptake.  

why we can't make the kl calculated on the fly? 
monthly value probably the best you could have but suspect that it still won't reflect the real world senarios. again, temperature and open canopy.  

Slurp is probably the best way to go so far, since it accounts the incomplete canopy and tempareture. 
curving fitting is not recommended since the assumptions are very easy to be violated in NZ.  

#### 20200330 Grab above ground component and soil water 

Unit conversion for visulising the graph in apsimx UI.   

```{r, include=FALSE}
con <- dbConnect(SQLite(), here::here("01Data/ProcessedData/Richard.sqlite3"))
dbListTables(con)
# Bring in the meta data for calculation 
meta_data <- dbReadTable(con, "meta_data") %>% 
  as.data.table()
biomass <- dbReadTable(con, "biomass") %>% 
  as.data.table()
water <- dbReadTable(con, "SoilWater", check.names = FALSE) %>% 
  as.data.table()
root <- dbReadTable(con, "root") %>% 
  as.data.table()
ES <- dbReadTable(con, "ES") %>% 
  as.data.table()
dbDisconnect(con)

# Fix the date 
biomass[, Clock.Today := as.Date(Clock.Today, origin = "1970-01-01")]
water[, Clock.Today := as.Date(Clock.Today, origin = "1970-01-01")]



# Dry matter value to evaluate the apsimx predict value 
# Need to be in g/m2 unit
DM_obs <- meta_data %>% 
  dplyr::filter(!is.na(Equations)) %>% 
  dplyr::mutate(Equations = gsub("\\s|\\(|\\)",".", Equations))
```

## Evapotranspiration surface

This is for calibrate the model for Potential Ep since there is no crop growing on this area?


```{r ES, fig.height=7, fig.width=7, dpi=300, echo=FALSE, message=FALSE, warning=FALSE, fig.cap="Surface evaporation data."}

vars <- colnames(ES[, Data:= NULL])[c(1:6)]
ES_melted <- ES %>% 
  melt(.,id.vars = vars, variable.factor = FALSE) %>% 
  .[, variable := as.integer(gsub("SWC.0\\.", "", variable))]
ES %>% 
  ggplot(aes(Clock.Today, SWCES_0.0.2)) + 
  geom_point() +
  facet_grid(  Experiment ~ .) +
  theme_water()

```

## Biomass Overview 

**The Biomass has `r dim(biomass)[1]` obs for two sites. **

```{r biomass, fig.height=12, fig.width=12, message=FALSE,warning=FALSE,dpi=300}

nums <- inspect_num(biomass)
show_plot(nums)
inspect_cat(biomass) %>% 
  show_plot()
```

```{r}
# GGally::ggpairs(biomass[,.(Site, Sowing.Date, Date, FW..g.,DW..g.)],aes(colour = Site, alpha = 0.4))
# colnames(biomass)
biomass_sub <- biomass[,.(Experiment, Clock.Today, SowingDate, DAS,  Shoot.No., Shoots.m2, FW..g., DW..g., Leaf.Area..cm.,LAImod)]
```

_Note:_

Too many vars, take long time to generate the pairplot. 
R is a bit slow to plot big correlation plots 
Try seaborn in python. still slow. 
consider to set up `Draker` plan. 

## Root Overview 


**The Root tab has `r dim(root)[1]` obs. **


```{r root, fig.height=12, fig.width=12,dpi=300}

nums <- inspect_num(root)
show_plot(nums)
```