---
title: "7. Point Pattern to Density"
author: "Wolfgang Hamer, Daniel Knitter"
output: html_document
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  warning = FALSE,
  message = FALSE,
  eval = TRUE,
  results = "show",
  #echo = FALSE,
  comment = "#>")
```


```{r libraries,warning=FALSE,message=FALSE}
library(magrittr)
library(knitr)
library(kableExtra)
library(sf)
library(raster)
library(FuzzyLandscapes)
library(dplyr)
library(rpredmod)
library(ggplot2)
library(spatstat)
library(maptools)
```

## Read point patterns as generated by "6_Fuzzy_to_Point_Pattern.Rmd"

```{r}
binary_addition_p <- readRDS(file = "../data/derived_data/vector/binary_addition_p.rds")
Fc_binary_addition <- readRDS(file = "../data/derived_data/vector/Fc_binary_addition.rds")
F_reindeer <- readRDS(file = "../data/derived_data/vector/F_reindeer.rds")
Fc_reindeer <- readRDS(file = "../data/derived_data/vector/Fc_reindeer.rds")
F_elk <- readRDS(file = "../data/derived_data/vector/F_elk.rds")
Fc_elk <- readRDS(file = "../data/derived_data/vector/Fc_elk.rds")
Findings <- readRDS(file = "../data/derived_data/vector/Findings.rds")
Findings_reindeer <- readRDS(file = "../data/derived_data/vector/Findings_reindeer.rds")
Findings_elk <- readRDS(file = "../data/derived_data/vector/Findings_elk.rds")
```


## Calculate the optimal sigma value using function `bw.ppl`

```{r}
(opt_sigma_vals <- data.frame(variable = c("F_reindeer",
                                           "Fc_reindeer",
                                           "F_elk",
                                           "Fc_elk",
                                           "Findings",
                                           "Findings_reindeer",
                                           "Findings_elk"),
                              opt_sig = c(bw.ppl(F_reindeer),
                                          bw.ppl(Fc_reindeer),
                                          bw.ppl(F_elk),
                                          bw.ppl(Fc_elk),
                                          bw.ppl(Findings),
                                          bw.ppl(Findings_reindeer),
                                          bw.ppl(Findings_elk))))

(sigmaval <- median(opt_sigma_vals$opt_sig))
resolu <- 250
```

## Create density raster

```{r}
reindeer_rand_dens <- density(F_reindeer, sigma=sigmaval,eps=c(resolu,resolu))
reindeer_clust_dens <- density(Fc_reindeer, sigma=sigmaval,eps=c(resolu,resolu))

elk_rand_dens <- density(F_elk, sigma=sigmaval,eps=c(resolu,resolu))
elk_clust_dens <- density(Fc_elk, sigma=sigmaval,eps=c(resolu,resolu))

Findings_dens <- density(Findings, sigma=sigmaval,eps=c(resolu,resolu))
Findings_reindeer_dens <- density(Findings_reindeer, sigma=sigmaval,eps=c(resolu,resolu))
Findings_elk_dens <- density(Findings_elk, sigma=sigmaval,eps=c(resolu,resolu))
```

```{r}
dens_stack <- stack(raster(reindeer_rand_dens), 
                    raster(reindeer_clust_dens),
                    raster(elk_rand_dens), 
                    raster(elk_clust_dens),
                    raster(Findings_reindeer_dens),
                    raster(Findings_elk_dens))
names(dens_stack) <- c("reindeer_rand_dens", 
                       "reindeer_clust_dens",
                       "elk_rand_dens", 
                       "elk_clust_dens",
                       "Findings_reindeer_dens",
                       "Findings_elk_dens")
crs(dens_stack) <- CRS("+init=epsg:25832")

plot(dens_stack,nc=2)

writeRaster(dens_stack, 
            filename=paste0("../data/derived_data/raster/density/raw/",names(dens_stack)), 
            bylayer=TRUE,
            format="GTiff",
            overwrite=TRUE)
```

## Crop to study area

```{r, fig.width=10}
sh <- st_read("../data/derived_data/vector/sh_nomarsh.shp")
cropstack <- mask(dens_stack, sh)
```

## Rescale data by range 

```{r}
scaled_cropstack <- stack(
  Map(function(lay){
    (1/diff(range(lay[],na.rm=TRUE)))*(lay-min(lay[],na.rm=TRUE))
  },
  lay=as.list(cropstack)))

writeRaster(scaled_cropstack, 
            filename=paste0("../data/derived_data/raster/density/scaled/",names(scaled_cropstack)),
            bylayer=TRUE,
            format="GTiff",
            overwrite=TRUE)
```

## Calculate the difference of the fuzzy densities to Findings density

```{r, fig.height= 3}
plot(stack(scaled_cropstack$Findings_reindeer_dens,scaled_cropstack$Findings_elk_dens))

diffstack <- stack(
  (scaled_cropstack[[which(grepl("rein",names(scaled_cropstack)))]]-scaled_cropstack$Findings_reindeer_dens) %>% 
    `names<-`(names(scaled_cropstack)[which(grepl("rein",names(scaled_cropstack)))]) %>% 
    dropLayer(.,3),
  (scaled_cropstack[[which(grepl("elk",names(scaled_cropstack)))]]-scaled_cropstack$Findings_elk_dens) %>% 
    `names<-`(names(scaled_cropstack)[which(grepl("elk",names(scaled_cropstack)))]) %>% 
    dropLayer(.,3))

writeRaster(diffstack, 
            filename=paste0("../data/derived_data/raster/density/difference/",names(diffstack)),
            bylayer=TRUE,
            format="GTiff",
            overwrite=TRUE)
```

```{r, fig.height= 7}
breakss <- c(-1,
             -0.75,
             -0.25,
             0.25,
             0.75,
             1)
colss <- c("#e66101","#fdb863","#1a9641","#b2abd2","#5e3c99")

plot(diffstack,
     nc=2,
     breaks = breakss,
     col = colss,
     legend=T,
     axis.args=list(at=breakss,
                    labels=paste(round(breakss,3)), 
                    cex.axis=.8))
```


