report_prototype.rmd

---
title: "EMMA Report Prototype"
description: Modeling vegetation postfire recovery data
author:
  - name: Adam Wilson & Glenn Moncrieff
editor_options: 
  chunk_output_type: console
output:
    html_document:
      toc: true
      toc_depth: 2
---

Model last updated at `r now()`.


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

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

# Model Overview
The details are given in [@slingsby_near-real_2020;@wilson_climatic_2015], but in short what we do is 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. For this we use a negative exponential curve with the following form:

$$\mu_{i,t}=\alpha_i+\gamma_i\Big(1-e^{-\frac{age_{i,t}}{\lambda_i}}\Big)$$

where $\mu_{i,t}$ is the expected NDVI for site $i$ at time $t$

The observed greenness $NDVI_{i,t}$  is assumed to follow a normal distribution with mean $\mu_{i,t}$
$$NDVI_{i,t}\sim\mathcal{N}(\mu_{i,t},\sigma_)$$

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. The full model also includes a sinusoidal term to capture seasonal variation, but lets keep it simple here.  


## Workflow

This repository was developed using the Targets framework as follows.


```{r make_workfow, echo=F, eval = TRUE, message = F, include=FALSE}

#This code chunk makes the targets visualization but does not call it (see next chunk for that)

#tfile=paste0(tempfile(),".html")
  targets_workflow_plot  <-targets::tar_visnetwork(reporter = "silent",
                            targets_only = TRUE,) #%>% 
#  htmlwidgets::saveWidget(file = tfile)
#webshot::install_phantomjs()
#webshot::webshot(tfile, "network.png")
#![targets_network](network.png)

```

```{r plot_workflow, echo=F, eval = TRUE, message = F, fig.width = 10, fig.height = 6}

  #This code chunk plots the targets visualization.
  # I do this in two steps to avoid errors that otherwise can occur

  targets_workflow_plot

```

## Results

### Environmental Controls on Ecosystem Recovery

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 = TRUE, 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_line(aes(y=median),colour="blue") +
  geom_line(aes(y=y_obs),colour="black",lwd=0.5,alpha=0.3) +
  geom_ribbon(aes(ymin=q5,ymax=q95),alpha=0.5)+
  facet_wrap(~cellID) +
  labs(x="time since fire (years)",y="NDVI") +
  theme_bw()

```


## Spatial Predictions

Maps of spatial parameters in the model.  

```{r compare_data2, echo = F, eval = TRUE, warning = F, message = FALSE}

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")

leaflet() %>% setView(lng = 18.577, lat = -33.998707, zoom = 10) %>% 
  addProviderTiles(providers$Esri.WorldImagery) %>% 
  addRasterImage(rast[[1]],group=names(rast)[1]) %>% #, color = ~pal(values(rast)[,1]) 
  addRasterImage(rast[[2]],group=names(rast)[2]) %>% 
  addRasterImage(rast[[3]],group=names(rast)[3]) %>% 
  addRasterImage(rast[[4]],group=names(rast)[4]) %>% 
  addLayersControl(
    baseGroups = names(rast),
    options = layersControlOptions(collapsed = FALSE)) #%>% 
#  addLegend("bottomright", pal = pal, values = ~,
#    title = "Est. GDP (2010)",
#    labFormat = labelFormat(prefix = "$"),
#    opacity = 1
#  )
```


# Park-specific Information

### NDVI and Time Since Fire

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

#CREATE wgs84 version of the park polygon

    focal_park %>% 
    st_transform(crs = st_crs(4326)) -> focal_wgs



library(leaflet)

#Note that the "addLegend_decreasing" isn't quite right, but might be close enough if you feel strongly about the ordering

addLegend_decreasing <- function (map, position = c("topright", "bottomright", "bottomleft", 
			    "topleft"), pal, values, na.label = "NA", bins = 7, colors, 
		  opacity = 0.5, labels = NULL, labFormat = labelFormat(), 
		  title = NULL, className = "info legend", layerId = NULL, 
		  group = NULL, data = getMapData(map), decreasing = FALSE) {
	position <- match.arg(position)
	type <- "unknown"
	na.color <- NULL
	extra <- NULL
	if (!missing(pal)) {
		if (!missing(colors)) 
			stop("You must provide either 'pal' or 'colors' (not both)")
		if (missing(title) && inherits(values, "formula")) 
			title <- deparse(values[[2]])
		values <- evalFormula(values, data)
		type <- attr(pal, "colorType", exact = TRUE)
		args <- attr(pal, "colorArgs", exact = TRUE)
		na.color <- args$na.color
		if (!is.null(na.color) && col2rgb(na.color, alpha = TRUE)[[4]] == 
		    0) {
			na.color <- NULL
		}
		if (type != "numeric" && !missing(bins)) 
			warning("'bins' is ignored because the palette type is not numeric")
		if (type == "numeric") {
			cuts <- if (length(bins) == 1) 
				pretty(values, bins)
			else bins	
			
			if (length(bins) > 2) 
				if (!all(abs(diff(bins, differences = 2)) <= 
				         sqrt(.Machine$double.eps))) 
					stop("The vector of breaks 'bins' must be equally spaced")
			n <- length(cuts)
			r <- range(values, na.rm = TRUE)
			cuts <- cuts[cuts >= r[1] & cuts <= r[2]]
			n <- length(cuts)
			p <- (cuts - r[1])/(r[2] - r[1])
			extra <- list(p_1 = p[1], p_n = p[n])
			p <- c("", paste0(100 * p, "%"), "")
			if (decreasing == TRUE){
				colors <- pal(rev(c(r[1], cuts, r[2])))
				labels <- rev(labFormat(type = "numeric", cuts))
			}else{
				colors <- pal(c(r[1], cuts, r[2]))
				labels <- rev(labFormat(type = "numeric", cuts))
			}
			colors <- paste(colors, p, sep = " ", collapse = ", ")
			
		}
		else if (type == "bin") {
			cuts <- args$bins
			n <- length(cuts)
			mids <- (cuts[-1] + cuts[-n])/2
			if (decreasing == TRUE){
				colors <- pal(rev(mids))
				labels <- rev(labFormat(type = "bin", cuts))
			}else{
				colors <- pal(mids)
				labels <- labFormat(type = "bin", cuts)
			}
			
		}
		else if (type == "quantile") {
			p <- args$probs
			n <- length(p)
			cuts <- quantile(values, probs = p, na.rm = TRUE)
			mids <- quantile(values, probs = (p[-1] + p[-n])/2, 
				 na.rm = TRUE)
			if (decreasing == TRUE){
				colors <- pal(rev(mids))
				labels <- rev(labFormat(type = "quantile", cuts, p))
			}else{
				colors <- pal(mids)
				labels <- labFormat(type = "quantile", cuts, p)
			}
		}
		else if (type == "factor") {
			v <- sort(unique(na.omit(values)))
			colors <- pal(v)
			labels <- labFormat(type = "factor", v)
			if (decreasing == TRUE){
				colors <- pal(rev(v))
				labels <- rev(labFormat(type = "factor", v))
			}else{
				colors <- pal(v)
				labels <- labFormat(type = "factor", v)
			}
		}
		else stop("Palette function not supported")
		if (!any(is.na(values))) 
			na.color <- NULL
	}
	else {
		if (length(colors) != length(labels)) 
			stop("'colors' and 'labels' must be of the same length")
	}
	legend <- list(colors = I(unname(colors)), labels = I(unname(labels)), 
	               na_color = na.color, na_label = na.label, opacity = opacity, 
	               position = position, type = type, title = title, extra = extra, 
	               layerId = layerId, className = className, group = group)
	invokeMethod(map, data, "addLegend", legend)
}


# Create continuous palettes

  pal <- colorNumeric(palette = "Reds",
                      domain = fires_wgs$Years,
                      reverse = TRUE)

  antipal <- colorNumeric(palette = "Reds",
                          domain = fires_wgs$Years,reverse = FALSE)
  
  
  ndvi_pal <- colorNumeric(palette = c( "#FFFFCC", "#41B6C4","#0C2C84"),
                           domain =  values(most_recent_ndvi_raster),
                           na.color = "transparent")

  

# Create bounding box for plotting  
  
bbox <- st_bbox(focal_wgs) %>%
  as.vector()



```


```{r ndvi_and_burns, fig.width = 10, fig.height = 7, echo = FALSE, warning = FALSE, message = FALSE}

  leaflet(data = focal_wgs) %>%
  addProviderTiles("Esri.NatGeoWorldMap", group = "NatGeo") %>%
    #addProviderTiles("NASAGIBS.ModisTerraTrueColorCR", group = "True Colors") %>%
    addProviderTiles(providers$Esri.WorldImagery, group = "World Imagery") %>%
  addPolygons(color = "black",
              stroke = TRUE,
              fill = FALSE,
              group = "Park") %>%
  addRasterImage(x = raster::raster(most_recent_ndvi_raster),
                 group = "NDVI")%>%
  addPolygons(data = fires_wgs,
              color = "red",
              fill = TRUE,fillOpacity = 0,
              stroke = TRUE,
              label = paste(fires_wgs$Years, "years since burn"),
              weight = 1,
              group = "Fire Outlines") %>%
    addPolygons(data = fires_wgs,
              color = ~antipal(Years),
              fillOpacity = .75,
              stroke = FALSE,
              label = fires_wgs$Years,
              group = "Fires") %>%
    addLegend_decreasing(position = "bottomright",
            pal = pal,
            values = ~fires_wgs$Years,
    title = "Years<br/>Since<br/>Fire",
    opacity = 1)%>%
  addLegend(position = "bottomright",
            pal = ndvi_pal,          
            values = ~values(raster::raster(most_recent_ndvi_raster)),
            opacity = 1,
            title = "NDVI") %>%
    addLayersControl(
    baseGroups = c("NatGeo","World Imagery"),
    overlayGroups = c("NDVI", "Fire Outlines","Fires", "Park"),
    options = layersControlOptions(collapsed = FALSE),position = "topright") %>%
  fitBounds(bbox[1], bbox[2], bbox[3], bbox[4])

```

Figure X. Leaflet map showing NDVI and the number of years since the most recent fire.  NA values in the NDVI layer indicate missing data (e.g. due to clouds, water, etc.).  NA values in the fire layer indicate that no fires have been recorded for that area.



```{r basemap1, eval=FALSE, fig.width = 10, echo = FALSE, warning = FALSE, message = FALSE}

# Time since fire
  #library(ggmmap) # doesn't play well with ggplot
  #library(basemaps) # doesn't play well with ggplot

  #Trying with basemapR
  #library(basemapR)


  suppressMessages(
  ggplot()+
    base_map(bbox = st_bbox(focal_wgs), basemap = "google-terrain",nolabels = TRUE,increase_zoom = 2)+
    geom_sf(data = focal_wgs, fill=NA)+
    geom_sf(data = fires_wgs,
            mapping = aes(fill = Years),alpha=0.5)+
    coord_sf(xlim = st_bbox(focal_wgs)[c(1,3)],
           ylim = st_bbox(focal_wgs)[c(2,4)])+
    scale_fill_gradient(low="white", high="red",name = "Years \nSince \nBurn")
  )
  
  

```

```{r basemap2, eval=FALSE, fig.width = 10, echo = FALSE, warning = FALSE, message = FALSE}

# Time since fire
  #library(ggmmap) # doesn't play well with ggplot
  #library(basemaps) # doesn't play well with ggplot

  #Trying with basemapR
  ##devtools::install_github('Chrisjb/basemapR')
  #library(basemapR)

    focal_park %>% 
    st_transform(crs = st_crs(4326)) -> focal_wgs

  suppressMessages(
  ggplot()+
    base_map(bbox = st_bbox(focal_wgs), basemap = "google-satellite",increase_zoom = 2)+
    geom_sf(data = focal_wgs, fill=NA)+
    geom_sf(data = fires_wgs,
            mapping = aes(fill = Years),alpha=0.5)+
    coord_sf(xlim = st_bbox(focal_wgs)[c(1,3)],
           ylim = st_bbox(focal_wgs)[c(2,4)])+
    scale_fill_gradient(low="white", high="red",name = "Years \nSince \nBurn")
  )
  
  
  
```


### Climate Data

```{r precip, fig.width = 10, echo = FALSE, warning = FALSE, message = FALSE}

#construct needed bounding box
  bb1 <- focal_wgs %>%st_bbox()
  class(bb1) <- "numeric"

# get climate data.  Currently set to grab climate data starting from one day before, using 2 days if that fails.
  #If we fun into errors in the future, this could be turned into a "while" loop.
  
  
focal_arc <- 
  tryCatch(expr = arc2(date = seq(Sys.Date() - (time_window_days+1), Sys.Date() - 1, "days"),
                    box = bb1) %>%
                bind_rows(), error = function(e){
                  
                  arc2(date = seq(Sys.Date() - (time_window_days+2), Sys.Date() - 2, "days"),
                    box = bb1) %>%
                bind_rows()

                })


#in case the park is too small to show up, enlarge the bounding box until some data shows up  

  if(nrow(focal_arc) == 0){
    
    while(nrow(focal_arc)==0){
      
      bb1[1] <- bb1[1] - 0.5
      bb1[2] <- bb1[2] - 0.5
      bb1[3] <- bb1[3] + 0.5
      bb1[4] <- bb1[4] + 0.5
      
      
      focal_arc <- 
  tryCatch(expr = arc2(date = seq(Sys.Date() - (time_window_days+1), Sys.Date() - 1, "days"),
                    box = bb1) %>%
                bind_rows(), error = function(e){
                  
                  arc2(date = seq(Sys.Date() - (time_window_days+2), Sys.Date() - 2, "days"),
                    box = bb1) %>%
                bind_rows()

                })
      
    }
    
    
    
  }#if no data retereived





  
#toss values outside of the park (unless there aren't any)  

  
  focal_arc %>%
    st_as_sf(coords = c("lon","lat"),
             crs = st_crs(focal_wgs)) -> focal_arc

st_intersection(x = focal_arc,
                y = focal_wgs) -> park_arc


#make a plot
  if(nrow(park_arc) == 0){
    
    
    
    focal_arc %>%
      group_by(date)%>%
      summarize(mean_precip = mean(precip)) %>%
      st_drop_geometry()%>% mutate(date = as.Date(date))%>%
      ggplot(mapping = aes(x = date, y = mean_precip))+
      geom_line()+
      xlab("Date")+
      ylab("Mean Precipitation (mm)")+
             ylim(c(0, max(focal_arc$precip)%>% {if(.==0){100}else .}) )
      

  }else{
    
    park_arc %>%
      group_by(date)%>%
      summarize(mean_precip = mean(precip)) %>%
      st_drop_geometry()%>% mutate(date = as.Date(date))%>%
      ggplot(mapping = aes(x = date, y = mean_precip))+
      geom_line()+
      xlab("Date")+
      ylab("Mean Precipitation (mm)")+
             ylim(c(0, max(focal_arc$precip)%>% {if(.==0){100}else .}) )
  
  }



  #clean up
    rm(focal_arc, park_arc)

  #clear out arc data
    arc2_cache$delete_all()

```