26-chillR_eval_PLS.Rmd

# Evaluating PLS outputs  {#pls_eval_PLS}

## Learning goals for this lesson {-#goals_PLS_eval}

- Understand what we can do with the temperature response phase delineations
- Be able to compute estimates of chill and heat requirements based on the outputs of a PLS analysis
- Be able to make temperature response plots for phenology data

## Chilling and forcing requirements

Despite some trouble in delineating chill and heat response periods, we've often been able to identify both dormancy phases. Let's see how we can make use of this information. We'll use the example of our 'Alexander Lucas' pears again, for which we'll first draw the PLS response curve.

```{r, echo=FALSE}
library(kableExtra)
```


```{r}
library(chillR)
Alex_first <- read_tab("data/Alexander_Lucas_bloom_1958_2019.csv") %>%
  select(Pheno_year, First_bloom) %>%
  mutate(Year = as.numeric(substr(First_bloom, 1, 4)),
         Month = as.numeric(substr(First_bloom, 5, 6)),
         Day = as.numeric(substr(First_bloom, 7, 8))) %>%
  make_JDay() %>%
  select(Pheno_year, 
         JDay) %>%
  rename(Year = Pheno_year,
         pheno = JDay)

temps <- read_tab("data/TMaxTMin1958-2019_patched.csv")

temps_hourly <- temps %>%
  stack_hourly_temps(latitude = 50.6)

daychill <- daily_chill(hourtemps = temps_hourly,
                        running_mean = 1,
                        models = list(Chilling_Hours = Chilling_Hours,
                                      Utah_Chill_Units = Utah_Model,
                                      Chill_Portions = Dynamic_Model,
                                      GDH = GDH)
                        )


plscf <- PLS_chill_force(daily_chill_obj = daychill,
                         bio_data_frame = Alex_first,
                         split_month = 6,
                         chill_models = "Chill_Portions",
                         heat_models = "GDH",
                         runn_means = 11)

```

We'll plot the results using the `plot_PLS_chill_force` function we produced in [PLS regression with agroclimatic metrics].


```{r, echo=FALSE}
library(ggplot2)
plot_PLS_chill_force <- function(plscf,
                                 chill_metric = "Chill_Portions",
                                 heat_metric = "GDH",
                                 chill_label = "CP",
                                 heat_label = "GDH",
                                 chill_phase = c(-48, 62),
                                 heat_phase = c(3, 105.5))
{
  PLS_gg <- plscf[[chill_metric]][[heat_metric]]$PLS_summary %>%
    mutate(Month = trunc(Date/100),
           Day = Date - Month * 100,
           Date = ISOdate(2002,
                          Month,
                          Day))
  
  PLS_gg[PLS_gg$JDay <= 0,"Date"]<-
    ISOdate(2001,
            PLS_gg$Month[PLS_gg$JDay <= 0],
            PLS_gg$Day[PLS_gg$JDay <= 0])
  
  PLS_gg <- PLS_gg %>%
    mutate(VIP_importance = VIP >= 0.8,
           VIP_Coeff = factor(sign(Coef) * VIP_importance))
  
  chill_start_date <- ISOdate(2001,
                              12,
                              31) + chill_phase[1] * 24 * 3600
  chill_end_date <- ISOdate(2001,
                            12,
                            31) + chill_phase[2] * 24 * 3600
  heat_start_date <- ISOdate(2001,
                             12,
                             31) + heat_phase[1] * 24 * 3600
  heat_end_date <- ISOdate(2001,
                           12,
                           31) + heat_phase[2] * 24 * 3600




  temp_plot <- ggplot(PLS_gg) +
    annotate("rect",
             xmin = chill_start_date,
             xmax = chill_end_date,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "blue") +
    annotate("rect",
             xmin = heat_start_date,
             xmax = heat_end_date,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "red") +
    annotate("rect",
             xmin = ISOdate(2001,
                            12,
                            31) +
               min(plscf$pheno$pheno,
                   na.rm = TRUE) * 24 * 3600,
             xmax = ISOdate(2001,
                            12,
                            31) +
               max(plscf$pheno$pheno,
                   na.rm = TRUE) * 24 * 3600,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "black") +
    geom_vline(xintercept = ISOdate(2001,
                                    12,
                                    31) +
                 median(plscf$pheno$pheno,
                        na.rm=TRUE) * 24 * 3600,
               linetype = "dashed") +
    geom_ribbon(aes(x = Date,
                    ymin = MetricMean - MetricStdev ,
                    ymax = MetricMean + MetricStdev ),
                fill = "grey") +
    geom_ribbon(aes(x = Date,
                    ymin = MetricMean - MetricStdev * (VIP_Coeff == -1),
                    ymax = MetricMean + MetricStdev * (VIP_Coeff == -1)),
                fill = "red") +
    geom_ribbon(aes(x = Date,
                    ymin = MetricMean - MetricStdev * (VIP_Coeff == 1),
                    ymax = MetricMean + MetricStdev * (VIP_Coeff == 1)),
                fill = "dark green") +
    geom_line(aes(x = Date,
                  y = MetricMean)) +
    facet_wrap(vars(Type),
               scales = "free_y",
               strip.position = "left",
               labeller = 
                 labeller(
                   Type =
                     as_labeller(c(Chill = paste0("Chill (",
                                                  chill_label,
                                                  ")"),
                                   Heat = paste0("Heat (",
                                                 heat_label,
                                                 ")"))))) +
    ggtitle("Daily chill and heat accumulation rates") +
    theme_bw(base_size = 15) + 
    theme(strip.background = element_blank(),
          strip.placement = "outside",
          strip.text.y = element_text(size = 12),
          plot.title = element_text(hjust = 0.5),
          axis.title.y = element_blank()
          )
  
  VIP_plot <- ggplot(PLS_gg,
                     aes(x = Date,
                         y = VIP)) +
    annotate("rect",
             xmin = chill_start_date,
             xmax = chill_end_date,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "blue") +
    annotate("rect",
             xmin = heat_start_date,
             xmax = heat_end_date,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "red") +
    annotate("rect",
             xmin = ISOdate(2001,
                            12,
                            31) + min(plscf$pheno$pheno,
                                      na.rm = TRUE) * 24 * 3600,
             xmax = ISOdate(2001,
                            12,
                            31) + max(plscf$pheno$pheno,
                                      na.rm = TRUE) * 24 * 3600,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "black") +
    geom_vline(xintercept = ISOdate(2001,
                                    12,
                                    31) + median(plscf$pheno$pheno,
                                                 na.rm = TRUE) * 24 * 3600,
               linetype = "dashed") +
    geom_bar(stat = 'identity',
             aes(fill = VIP>0.8)) +
    facet_wrap(vars(Type), 
               scales = "free",
               strip.position = "left",
               labeller = 
                 labeller(
                   Type = as_labeller(c(Chill="VIP for chill",
                                        Heat="VIP for heat")))) +
    scale_y_continuous(
      limits = c(0,
                 max(plscf[[chill_metric]][[heat_metric]]$PLS_summary$VIP))) +
    ggtitle("Variable Importance in the Projection (VIP) scores") +
    theme_bw(base_size = 15) + 
    theme(strip.background = element_blank(),
          strip.placement = "outside",
          strip.text.y = element_text(size = 12),
          plot.title = element_text(hjust = 0.5),
          axis.title.y = element_blank()
          ) +
    scale_fill_manual(name = "VIP", 
                      labels = c("<0.8", ">0.8"), 
                      values = c("FALSE" = "grey",
                                 "TRUE" = "blue")) +
    theme(axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.title.x = element_blank(),
          axis.title.y = element_blank())
  
  coeff_plot <- ggplot(PLS_gg,
                       aes(x = Date,
                           y = Coef)) +
    annotate("rect",
             xmin = chill_start_date,
             xmax = chill_end_date,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "blue") +
    annotate("rect",
             xmin = heat_start_date,
             xmax = heat_end_date,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "red") +
    annotate("rect",
             xmin = ISOdate(2001,
                            12,
                            31) + min(plscf$pheno$pheno,
                                      na.rm = TRUE) * 24 * 3600,
             xmax = ISOdate(2001,
                            12,
                            31) + max(plscf$pheno$pheno,
                                      na.rm = TRUE) * 24 * 3600,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "black") +
    geom_vline(xintercept = ISOdate(2001,
                                    12,
                                    31) + median(plscf$pheno$pheno,
                                                 na.rm = TRUE) * 24 * 3600,
               linetype = "dashed") +
    geom_bar(stat = 'identity',
             aes(fill = VIP_Coeff)) +
    facet_wrap(vars(Type),
               scales = "free",
               strip.position = "left",
               labeller =
                 labeller(
                   Type = as_labeller(c(Chill = "MC for chill",
                                        Heat = "MC for heat")))) +
    scale_y_continuous(
      limits = c(min(plscf[[chill_metric]][[heat_metric]]$PLS_summary$Coef),
                 max(plscf[[chill_metric]][[heat_metric]]$PLS_summary$Coef))) +
    ggtitle("Model coefficients (MC)") +
    theme_bw(base_size = 15) + 
    theme(strip.background = element_blank(),
          strip.placement = "outside",
          strip.text.y = element_text(size = 12),
          plot.title = element_text(hjust = 0.5),
          axis.title.y = element_blank()
          ) +
    scale_fill_manual(name = "Effect direction", 
                      labels = c("Advancing",
                                 "Unimportant",
                                 "Delaying"), 
                      values = c("-1" = "red",
                                 "0" = "grey",
                                 "1" = "dark green")) +
    ylab("PLS coefficient") +
    theme(axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.title.x = element_blank(),
          axis.title.y = element_blank())
  
  library(patchwork)
  
  plot <- (VIP_plot +
             coeff_plot +
             temp_plot +
             plot_layout(ncol = 1,
                         guides = "collect")
           ) & theme(legend.position = "right",
                     legend.text = element_text(size = 8),
                     legend.title = element_text(size = 10),
                     axis.title.x = element_blank())

plot

}



```


```{r}
plot_PLS_chill_force(plscf,
                     chill_metric = "Chill_Portions",
                     heat_metric = "GDH",
                     chill_label = "CP",
                     heat_label = "GDH",
                     chill_phase = c(-48, 62),
                     heat_phase = c(3, 105.5))


```

I've already indicated here what phases I consider relevant for chilling (light blue) and forcing (light red).

If we accept these periods as chilling and forcing phases, the amount of chill or heat accumulated during these phases should be a reasonable approximation of the tree's agroclimatic needs. We can easily calculate these with the `tempResponse` function. Of course, chill and heat accumulated during these phases isn't identical across all years of the dataset. It's not entirely clear how to arrive at reliable estimates for chill and heat requirements from this set of values, but I've found it quite reasonable to use the mean accumulation during the delineated periods as our best estimate and the standard deviation to evaluate the error.

```{r}

chill_phase <- c(317, 62)
heat_phase <- c(3, 105.5)

chill <- tempResponse(hourtemps = temps_hourly,
                      Start_JDay = chill_phase[1],
                      End_JDay = chill_phase[2],
                      models = list(Chill_Portions = Dynamic_Model),
                      misstolerance = 10)

heat <- tempResponse(hourtemps = temps_hourly,
                     Start_JDay = heat_phase[1],
                     End_JDay = heat_phase[2],
                     models = list(GDH = GDH))
```

Let's draw quick histograms for good measure and then quantify the chill and heat requirements, plus their errors.

```{r, warning=FALSE, message=FALSE}

ggplot(data = chill,
       aes(x = Chill_Portions)) +
  geom_histogram() +
  ggtitle("Chill accumulation during endodormancy (Chill Portions)") +
  xlab("Chill accumulation (Chill Portions)") +
  ylab("Frequency between 1958 and 2019") +
  theme_bw(base_size = 12)

ggplot(data = heat,
       aes(x = GDH)) +
  geom_histogram() +
  ggtitle("Heat accumulation during ecodormancy (GDH)") +
  xlab("Heat accumulation (Growing Degree Hours)") +
  ylab("Frequency between 1958 and 2019") +
  theme_bw(base_size = 12)


chill_requirement <- mean(chill$Chill_Portions)
chill_req_error <- sd(chill$Chill_Portions)

heat_requirement <- mean(heat$GDH)
heat_req_error <- sd(heat$GDH)

```

So we find that the chilling requirement (`chill_requirement`) is about `r round(chill_requirement,1)` Chill Portions, with an estimated error (`chill_req_error`) of `r round(chill_req_error,1)` CP. The heat requirement (`heat_requirement`) is about `r round(heat_requirement)` Growing Degree Hours, with an estimated error (`heat_req_error`) of `r round(heat_req_error)` GDH. Compared to other fruit trees, this is a fairly high chilling requirement, while the heat need is quite low.

I should add that we have some ambiguity here regarding the duration of the chilling period especially, and a considerable overlap between the two phases. This may lead to these estimates not being totally reliable here. There have been some discussions about how useful these statistical estimates are. In some cases, they have been pretty close to requirements that were estimated experimentally, while in other cases, there were substantial differences.

Even though some doubts about the accuracy of this method to determine agroclimatic needs remain, it has enabled quick and easy initial estimates of chill and heat requirements, at least for cultivars for which we have long-term bloom data (coupled with temperature recordings).

The first statistical estimates of chill and heat needs with this method were produced for 'Schneider's späte Knorpelkirsche' cherries in Klein-Altendorf [[@luedeling2013identification]](https://link.springer.com/article/10.1007/s00484-012-0594-y). In that study, we just used a temperature-based PLS analysis, but we later also ran similar analyses with PLS regression using agroclimatic metrics [[@luedeling2013differential]](https://www.sciencedirect.com/science/article/pii/S0168192313001780).

## Response to seasonal temperature

In addition to estimating chill and heat requirements, we can also take a more general look at tree responses to temperature during these delineated phases. While (as I hope you remember) our delineations of chilling and forcing phases were based on the assumption that the agroclimatic models were somewhat accurate, we'll now take a step back and use raw temperatures again. Now that we have a reasonable estimate of the major temperature response phases, let's see if we can explain (part of) the trees' bloom dates just by mean temperature during these periods.

We know that conditions during both the chilling and the forcing phases affect bloom dates, so let's try to illustrate both of these effects in one figure. We now have two independent variables (mean temperature during chilling and mean temperature during forcing), and we have one dependent variable (bloom date). This means we need three dimensions in our plot.

What `chillR` provides for us is a bit of a black box, but we'll reproduce this ourselves later. The function to look at for this is the `make_pheno_trend_plot` function. The only inputs we need to provide are a `data.frame` with weather data, a `data.frame` with phenology data and the start and end days of the temperature response phases we identified. The other inputs are necessary because, like many other functions in the package, `chillR` is currently set up to produce images rather than just plots in `R`.

```{r, eval=FALSE}

chill_phase <- c(317, 62)
heat_phase <- c(360, 106) # note that the end date here was rounded
                          # to an integer number, so that a proper 
                          # axis label can be generated.


mpt <- make_pheno_trend_plot(weather_data_frame = temps,
                             pheno = Alex_first,
                             Start_JDay_chill = chill_phase[1], 
                             End_JDay_chill = chill_phase[2],
                             Start_JDay_heat = heat_phase[1],
                             End_JDay_heat = heat_phase[2],
                             outpath = "data/",
                             file_name = "pheno_trend_plot",
                             plot_title =  "Impacts of chilling and forcing temperatures on pear phenology",
                             image_type = "png", 
                             colorscheme = "normal")

```


![Standard `chillR` output of phenology trend plot for pears 'Alexander Lucas' in Klein-Altendorf](pictures/pheno_trend_plot.png)

What we see here is a surface that is interpolated from all bloom dates included in the phenology record, which relates bloom date expectations to temperatures during the chilling phase (x-axis) and the forcing phase (y-axis). The shape of the plotted surface is defined by the array of winter conditions that have been observed since 1958. Apparently, for instance, there have been no years with an average temperature during chilling around -2°C and 4°C during forcing. Hence the plot area remains white in this region.

The array of historically observed conditions is relatively narrow here, with a fairly strong correlation between temperatures during eco- and endodormancy, as indicated by the narrow diagonal shape of the area. Nevertheless, we can already get an idea about the temperature responses. All the blue colors are near the top of the plot, while the red colors are at the bottom. This means that early bloom dates (blue) were mostly associated with a warm forcing phase, while cool conditions during that period were related to late bloom. Evaluating the effect of temperatures during chilling is more difficult here, because of missing data in climatic spaces that would be relevant. We'll look at other examples later that show a clearer response.

Now let's see if we can reproduce this figure using `ggplot`.

We'll first compute the mean temperatures for the chilling and forcing period. We've done something like this before, in the [Simple phenology analysis] chapter. I'll reuse code from there to make a function to summarize temperatures for a particular JDay interval

```{r}
mean_temp_period <- function(
    temps,
    start_JDay,
    end_JDay,
    end_season = end_JDay)
{ temps_JDay <- make_JDay(temps) %>%
  mutate(Season =Year)
  
  if(start_JDay > end_season)
    temps_JDay$Season[which(temps_JDay$JDay >= start_JDay)]<-
        temps_JDay$Year[which(temps_JDay$JDay >= start_JDay)]+1
  
  if(start_JDay > end_season)
    sub_temps <- subset(temps_JDay,
                        JDay <= end_JDay | JDay >= start_JDay)
  
  if(start_JDay <= end_JDay)
    sub_temps <- subset(temps_JDay,
                        JDay <= end_JDay & JDay >= start_JDay)
  
  mean_temps <- aggregate(sub_temps[, c("Tmin", "Tmax")],
                          by = list(sub_temps$Season),
                          FUN = function(x) mean(x,
                                                 na.rm=TRUE))
  
  mean_temps[, "n_days"] <- aggregate(sub_temps[, "Tmin"],
                                      by = list(sub_temps$Season),
                                      FUN = length)[,2]
  
  mean_temps[, "Tmean"] <- (mean_temps$Tmin + mean_temps$Tmax) / 2
  mean_temps <- mean_temps[, c(1, 4, 2, 3, 5)]
  colnames(mean_temps)[1] <- "End_year"
  
  return(mean_temps)
}

mean_temp_chill <- mean_temp_period(temps = temps,
                                    start_JDay = chill_phase[1],
                                    end_JDay = chill_phase[2],
                                    end_season = 60)

mean_temp_heat <- mean_temp_period(temps = temps,
                                   start_JDay = heat_phase[1],
                                   end_JDay = heat_phase[2],
                                   end_season = 60)

```

Note that we have to be a bit careful here that we don't assign heat or chill accumulation to the wrong season. This can easily happen when the chilling phase begins and ends before the beginning of January. To prevent this, I added the parameter `end_season` to the function. If we set this to a JDay of the year that we expect the phenological stage of interest to occur in, we should have no problem with erroneous assignment to seasons.

Now we can combine these two datasets (just the `Tmean` columns) and merge them with the phenology data. Before that, all incomplete years should be filtered out.

```{r}
mean_temp_chill <- 
  mean_temp_chill[which(mean_temp_chill$n_days >=
                          max(mean_temp_chill$n_days)-1),]
mean_temp_heat <-
  mean_temp_heat[which(mean_temp_heat$n_days >=
                         max(mean_temp_heat$n_days)-1),]

mean_chill <- mean_temp_chill[, c("End_year",
                                  "Tmean")]
colnames(mean_chill)[2] <- "Tmean_chill"

mean_heat <- mean_temp_heat[,c("End_year",
                               "Tmean")]
colnames(mean_heat)[2] <- "Tmean_heat"

phase_Tmeans <- merge(mean_chill,
                      mean_heat, 
                      by = "End_year")


pheno <- Alex_first
colnames(pheno)[1] <- "End_year"

Tmeans_pheno <- merge(phase_Tmeans,
                      pheno,
                      by = "End_year")

```

```{r, eval=FALSE}
head(Tmeans_pheno)
```


```{r,echo=FALSE}

kable(head(Tmeans_pheno)) %>%
  kable_styling("striped", position = "left", font_size=10)

```

Now we want to interpolate a surface from these data. There are various ways to do this. In spatial contexts, the so-called [Kriging](https://en.wikipedia.org/wiki/Kriging) technique is a widely used interpolation method. Check out the link for more details. `chillR` uses this technique as implemented in the [`fields`](https://cran.r-project.org/web/packages/fields/index.html) package (function `Krig`).

```{r, message=FALSE, warning=FALSE}

library(fields)
k <- Krig(x = as.matrix(
                Tmeans_pheno[,
                             c("Tmean_chill",
                               "Tmean_heat")]),
          Y = Tmeans_pheno$pheno)

pred <- predictSurface(k)
colnames(pred$z) <- pred$y
rownames(pred$z) <- pred$x

library(reshape2)
melted <- melt(pred$z)
  
library(metR)
library(colorRamps)
  
colnames(melted) <- c("Tmean_chill",
                      "Tmean_heat",
                      "value")


ggplot(melted,
       aes(x = Tmean_chill,
           y = Tmean_heat,
           z = value)) +
  geom_contour_fill(bins = 100) +
  scale_fill_gradientn(colours = alpha(matlab.like(15)),
                       name = "Bloom date \n(day of the year)") +
  geom_contour(col = "black")  +
  geom_point(data = Tmeans_pheno,
             aes(x = Tmean_chill,
                 y = Tmean_heat,
                 z = NULL),
             size = 0.7) +
  geom_text_contour(stroke = 0.2) +
  ylab(expression(paste("Forcing phase ", 
                        T[mean],
                        " (",
                        degree,
                        "C)"))) +
  xlab(expression(paste("Chilling phase ",
                        T[mean],
                        " (",
                        degree,
                        "C)")))  +
  theme_bw(base_size = 15)
    


```

This now looks remarkably easy, but it took me a while to put all the pieces together. What happened here is that we used the Kriging function and predicted a surface based on the result. The output contained separate elements for the x, y and z values of the surface, which had to be integrated into one `data.frame` to be useful for `ggplot2`. The `melt` function of the [`reshape2`](https://cran.r-project.org/web/packages/reshape2/index.html) package came in handy here. To make a proper plot from these data, we needed the [`colorRamps`](https://cran.r-project.org/web/packages/colorRamps/index.html) package for nice colors, and the [`metR`](https://cran.r-project.org/web/packages/metR/index.html) package to place nice labels on the contour lines.

Now that we have the code together to make nice plots, we should make a function from it.

```{r, message=FALSE,warning=FALSE}


pheno_trend_ggplot <- function(temps,
                               pheno,
                               chill_phase,
                               heat_phase,
                               phenology_stage = "Bloom")
{
  library(fields)
  library(reshape2)
  library(metR)
  library(ggplot2)
  library(colorRamps)
  
  # first, a sub-function (function defined within a function) to
  # compute the temperature means
  
  mean_temp_period <- function(temps,
                               start_JDay,
                               end_JDay, 
                               end_season = end_JDay)
    { temps_JDay <- make_JDay(temps) %>%
      mutate(Season = Year)

    if(start_JDay > end_season)
      temps_JDay$Season[which(temps_JDay$JDay >= start_JDay)] <-
        temps_JDay$Year[which(temps_JDay$JDay >= start_JDay)]+1
    
    if(start_JDay > end_season)
      sub_temps <- subset(temps_JDay,
                          JDay <= end_JDay | JDay >= start_JDay)
    
    if(start_JDay <= end_JDay)
      sub_temps <- subset(temps_JDay,
                          JDay <= end_JDay & JDay >= start_JDay)
    
    mean_temps <- aggregate(sub_temps[,
                                      c("Tmin",
                                        "Tmax")],
                            by = list(sub_temps$Season),
                            FUN = function(x) mean(x,
                                                   na.rm = TRUE))
    mean_temps[, "n_days"] <- aggregate(sub_temps[,
                                                  "Tmin"],
                                        by = list(sub_temps$Season),
                                        FUN = length)[,2]
    
    mean_temps[,"Tmean"] <- (mean_temps$Tmin + mean_temps$Tmax) / 2
    mean_temps <- mean_temps[, c(1, 4, 2, 3, 5)]
    colnames(mean_temps)[1] <- "End_year"
    return(mean_temps)
    }
  
  mean_temp_chill <- mean_temp_period(temps = temps,
                                      start_JDay = chill_phase[1],
                                      end_JDay = chill_phase[2],
                                      end_season = heat_phase[2])
  
  mean_temp_heat <- mean_temp_period(temps = temps,
                                     start_JDay = heat_phase[1],
                                     end_JDay = heat_phase[2],
                                     end_season = heat_phase[2])
  
  mean_temp_chill <-
    mean_temp_chill[which(mean_temp_chill$n_days >= 
                            max(mean_temp_chill$n_days)-1),]
  mean_temp_heat <-
    mean_temp_heat[which(mean_temp_heat$n_days >= 
                           max(mean_temp_heat$n_days)-1),]
  mean_chill <- mean_temp_chill[, c("End_year",
                                    "Tmean")]
  colnames(mean_chill)[2] <- "Tmean_chill"
  mean_heat<-mean_temp_heat[,c("End_year",
                               "Tmean")]
  colnames(mean_heat)[2] <- "Tmean_heat"
  phase_Tmeans <- merge(mean_chill,
                        mean_heat,
                        by = "End_year")
  
  colnames(pheno) <- c("End_year",
                       "pheno")
  Tmeans_pheno <- merge(phase_Tmeans,
                        pheno, 
                        by="End_year")
  
  # Kriging interpolation
  k <- Krig(x = as.matrix(Tmeans_pheno[,c("Tmean_chill",
                                          "Tmean_heat")]),
            Y = Tmeans_pheno$pheno)
  pred <- predictSurface(k)
  colnames(pred$z) <- pred$y
  rownames(pred$z) <- pred$x
  melted <- melt(pred$z)
  colnames(melted) <- c("Tmean_chill",
                        "Tmean_heat",
                        "value")
  
  ggplot(melted,
         aes(x = Tmean_chill,
             y = Tmean_heat,
             z = value)) +
    geom_contour_fill(bins = 60) +
    scale_fill_gradientn(colours = alpha(matlab.like(15)),
                         name = paste(phenology_stage,
                                      "date \n(day of the year)")) +
    geom_contour(col = "black") +
    geom_text_contour(stroke = 0.2) +
    geom_point(data = Tmeans_pheno,
               aes(x = Tmean_chill,
                   y = Tmean_heat,
                   z = NULL),
               size = 0.7)  +
    ylab(expression(paste("Forcing phase ",
                          T[mean],
                          " (",
                          degree,
                          "C)"))) +
    xlab(expression(paste("Chilling phase ",
                          T[mean],
                          " (",
                          degree,
                          "C)"))) +
    theme_bw(base_size = 15)
}


chill_phase <- c(317, 62)
heat_phase <- c(360, 105.5)

pheno_trend_ggplot(temps = temps,
                   pheno = Alex_first,
                   chill_phase = chill_phase,
                   heat_phase = heat_phase,
                   phenology_stage = "Bloom")

```

Now we've automated the whole procedure, so that we never have to bother with the details of the plot production again.

## Applying our functions to California walnuts

Let's use the functions we've produced so far to do a quick analysis of walnuts in California. With the following button, you can download leaf emergence data for the *Payne* walnut cultivar, observed in Davis, California. Save this in your `data` directory.

```{r, echo=FALSE, warning=FALSE}
Walnut_pheno <- read_tab("data/Davis_Payne_leaf_out.csv")
Walnut_pheno %>% download_this(
    output_name = "Davis_Payne_leaf_out",
    output_extension = ".csv",
    button_label = "Download walnut phenology data for Davis",
    button_type = "warning",
    has_icon = TRUE,
    icon = "fa fa-save"
  )
```

```{r, warning=FALSE, message=FALSE}
Cali_temps <- read_tab("data/Davis_weather.csv")
Walnut_pheno <- read_tab("data/Davis_Payne_leaf_out.csv") %>%
  mutate(Year = as.numeric(substr(Leaf.date,7,8)),
         Year = Year+(19+(Year<25))*100,
         Month = as.numeric(substr(Leaf.date,4,5)),
         Day = as.numeric(substr(Leaf.date,1,2))) %>%
  make_JDay() %>%
  select(Year, JDay)

colnames(Walnut_pheno) <- c("Year",
                            "pheno")

Cali_temps_hourly <- stack_hourly_temps(Cali_temps,
                                        latitude = 38.5)

Cali_daychill <- daily_chill(hourtemps = Cali_temps_hourly,
                             running_mean = 1,
                             models = list(Chilling_Hours = Chilling_Hours,
                                           Utah_Chill_Units = Utah_Model,
                                           Chill_Portions = Dynamic_Model,
                                           GDH = GDH)
    )


plscf <- PLS_chill_force(daily_chill_obj = Cali_daychill,
                         bio_data_frame = Walnut_pheno,
                         split_month = 6,
                         chill_models = "Chill_Portions",
                         heat_models = "GDH",
                         runn_means = 11)

plot_PLS_chill_force(plscf,
                     chill_metric = "Chill_Portions",
                     heat_metric = "GDH",
                     chill_label = "CP",
                     heat_label = "GDH",
                     chill_phase = c(-56, 5),
                     heat_phase = c(19, 77))

pheno_trend_ggplot(temps = Cali_temps,
                   pheno = Walnut_pheno,
                   chill_phase = c(309, 5),
                   heat_phase = c(19, 77),
                   phenology_stage = "Leaf emergence")

```

Here we have a more revealing color pattern. The earliest leaf emergence occurred when the chilling phase was cool, and the forcing phase warm. This is indicated by the blue colors in the top left corner of the plot. The bottom right corner, in contrast, looks red, indicating that a warm chilling phase and a cool forcing phase lead to late bloom.


## `Exercises` on evaluating PLS regression results {-#exercises_PLS_eval}

Please document all results of the following assignments in your `learning logbook`.

1) Reproduce the analysis for the 'Roter Boskoop' dataset.
2) We've looked at data from a number of locations so far. How would you expect this surface plot to look like in Beijing? And how should it look in Tunisia?