index.Rmd

---
title: "EMMA Prototype"
description: Modeling vegetation postfire recovery data
editor_options: 
  chunk_output_type: console
output:
  html_document:
    toc: true
    toc_depth: 2
---

```{r, echo=F, message=F,include=F, results="hide"}
library(targets)
library(tidyverse)
library(doParallel)
library(raster)
library(lubridate)
library(sf)
library(plotly)
library(leaflet)

# load data saved in the pipeline
tar_load(c(envdata, stan_data, model_results, spatial_outputs,model_prediction)) 
```

# Model Overview

We estimate the age of a site by calculating the years since the last fire. We then fit a curve to model the recovery of vegetation (measured using NDVI) as a function of it's age. An additional level models the parameters of the negative exponential curve as a function of environmental variables. This means that sites with similar environmental conditions should have similar recovery curves.

## Input data

```{r, echo=FALSE}

#get useful stats

tar_load(envdata_predict)
tar_load(stan_data_predict)

fitting_domain_pixels <- length(which(envdata$model_domain))

fitting_domain_long_pixels <- length(which(!is.na(envdata$long)& 
                                          envdata$model_domain))

predicting_domain_pixels <- length(which(envdata_predict$model_domain))

predicting_domain_long_pixels <- length(which(!is.na(envdata_predict$long)& 
                                          envdata_predict$model_domain))

rm(envdata_predict)

```



The model was last fit on `r now()`.

This version of the model was fit with `r filter(envdata,training) %>% nrow()` pixels including data from `r as_date(range(stan_data$y_date))[1]` to `r as_date(range(stan_data$y_date))[2]`. `r stan_data$J` out of `r predicting_domain_pixels` predicting domain pixels were used in model fitting (`r round(stan_data$J/predicting_domain_pixels,digits = 2)*100`%). Predictions were made for `r stan_data_predict$J` out of `r predicting_domain_pixels` prediction domain pixels (`r round(stan_data_predict$J/predicting_domain_pixels,digits = 2)*100`%).

## Workflow

This repository was developed using the Targets framework.

```{r, echo=F, eval=T, message=F, include=F, results="asis"}
#tfile=paste0(tempfile(),".html")
test<-targets::tar_visnetwork() #%>% 
#  htmlwidgets::saveWidget(file = tfile)
#webshot::install_phantomjs()
#webshot::webshot(tfile, "network.png")
#![targets_network](network.png)

```

```{r, echo=FALSE}
test
```

## Results

### Environmental Drivers

These parameters represent the relationship of the following environmental variables to the recovery trajectory.

```{r p1, echo=F, eval=T, warning=F, message=FALSE}

betas <- model_results %>%
  filter(type=="beta")


p1 <- ggplot(betas, aes(y=xname, xmin=q5,x=median,xmax=q95))+
  geom_pointrange(fill="grey")+
  facet_wrap(~parameter,nrow=1)+
  geom_vline(xintercept=0,col="grey")+
  xlab("Beta (regression coefficient +/- 95% CI)")+
  ylab("Environmental Variable") 

ggplotly(p1)

```

## Recovery Trajectories

The plot below illustrates some example recovery trajectories. It currently just shows the top 20 cells with the most observations.

```{r plot, echo=F, eval=T, message=FALSE,fig.height=12}

cells_with_long_records<-
  model_prediction %>% 
  group_by(cellID) %>% 
  summarize(n=n()) %>% 
  arrange(desc(n)) %>% 
  slice(1:20) # top 20 cells with the most observations

model_prediction %>% 
  filter(cellID%in%cells_with_long_records$cellID) %>% 
  ggplot(aes(x=age)) +
  geom_ribbon(aes(ymin=q5,ymax=q95),alpha=0.5)+
  geom_line(aes(y=y_obs),colour="darkred",lwd=0.5) +
  geom_line(aes(y=median),colour="blue") +
  facet_wrap(~cellID) +
  labs(x="time since fire (years)",y="NDVI") +
  theme_bw()

```

```{r, echo=F, eval=F}

cells_with_fires<-
  model_prediction %>% 
  filter(date>as_date("2010-01-01"),date<as_date("2015-01-01")) %>% 
  group_by(cellID) %>% 
  summarize(minage=min(age),n=n()) %>%
  filter(minage<1) %>% 
  arrange(desc(n)) %>% 
  slice(1:20) # top 20 cells with the most observations

model_prediction %>% 
  filter(cellID%in%cells_with_fires$cellID) %>% 
  ggplot(aes(x=date)) +
  geom_ribbon(aes(ymin=q5,ymax=q95),alpha=0.5)+
  geom_line(aes(y=y_obs),colour="darkred",lwd=0.5) +
  geom_line(aes(y=median),colour="blue") +
  facet_wrap(~cellID) +
  labs(x="time since fire (years)",y="NDVI") +
  theme_bw()

```

## Model Performance

Compare estimated vs observed values for all pixels. This is not true validation - these pixels were included in the model fitting.

```{r, message=F, echo=F}

model_prediction %>% 
  ggplot(aes(x=median,y=y_obs)) +
  geom_hex(bins=30)+
  geom_smooth(method = "lm",col="red")+
  geom_abline(color="blue")+
#  facet_wrap(~cellID) +
  scale_fill_viridis_c()+
  labs(x="Estimated NDVI",y="Observed NDVI",
       caption = "Blue line is 1:1, Red line is least squares regression. Count is the number of pixels in that location",
       title = "Estimated vs. Observed NDVI"
      ) +
  theme_bw()+
  coord_equal()

```


```{r, echo=FALSE}

tar_load(model_w_pred_summary_postfire_season_predict)

# Generating a plot of the alphas
model_w_pred_summary_postfire_season_predict %>%
  filter(grepl(pattern = "alpha[",fixed = TRUE,x=variable))-> alphas


alpha_rast <- spatial_outputs[[1]]
alpha_rast[1:ncell(alpha_rast)] <- NA
alpha_rast[stan_data$x_cellID] <- alphas$mean 
names(alpha_rast) <- "alpha" 

spatial_outputs <- stack(spatial_outputs,alpha_rast)

rm(model_w_pred_summary_postfire_season_predict)
```

## Spatial Predictions

Maps of spatial parameters in the model.

```{r, leaflet_map, echo=F, eval=T, warning=FALSE, message=FALSE,  results='asis'}

rast=projectRaster(spatial_outputs,crs = "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs +type=crs")

library(leafem)

leaflet() %>%
  setView(lng = 18.577, lat = -33.998707, zoom = 10) %>% 
  addProviderTiles(providers$Esri.WorldImagery) %>% 
  
  addRasterImage(rast[[1]],
                 group=names(rast)[1],
                 layerId = names(rast)[1]) %>% #, color = ~pal(values(rast)[,1]) 
  addImageQuery(rast[[1]],
                layerId = names(rast)[1]) %>% 
  
  addRasterImage(rast[[2]],
                 group=names(rast)[2],
                 layerId = names(rast)[2]) %>% 
    addImageQuery(rast[[2]],
                layerId = names(rast)[2]) %>% 

  addRasterImage(rast[[3]],group=names(rast)[3],
                 layerId = names(rast)[3]) %>%
    addImageQuery(rast[[3]],
                layerId = names(rast)[3]) %>% 

  addRasterImage(rast[[4]],group=names(rast)[4],
                 layerId = names(rast)[4]) %>%
    addImageQuery(rast[[4]],
                layerId = names(rast)[4]) %>% 

  addRasterImage(rast[[5]],group=names(rast)[5],
                 layerId = names(rast)[5]) %>%
    addImageQuery(rast[[5]],
                layerId = names(rast)[5]) %>% 

  addLayersControl(
    baseGroups = names(rast),
    options = layersControlOptions(collapsed = FALSE))  
#  addLegend("bottomright", pal = pal, values = ~,
#    title = "Est. GDP (2010)",
#    labFormat = labelFormat(prefix = "$"),
#    opacity = 1
#  )

```

# Park Reports

See the links below for park-level reports on vegetation status.

[Addo-Elephant_National_Park](reports/report.Addo-Elephant_National_Park.html)
[Agulhas_National_Park](reports/report.Agulhas_National_Park.html)
[Bontebok_National_Park](reports/report.Bontebok_National_Park.html)
[Garden_Route_National_Park](reports/report.Garden_Route_National_Park.html)
[Karoo_National_Park](reports/report.Karoo_National_Park.html)
[Namaqua_National_Park](reports/report.Namaqua_National_Park.html)
[Richtersveld_National_Park](reports/report.Richtersveld_National_Park.html)
[Table_Mountain_National_Park](reports/report.Table_Mountain_National_Park.html)
[Tankwa-Karoo_National_Park](reports/report.Tankwa-Karoo_National_Park.html)
[West_Coast_National_Park](reports/report.West_Coast_National_Park.html)