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”

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

(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))))
#>            variable   opt_sig
#> 1        F_reindeer 10951.019
#> 2       Fc_reindeer   967.041
#> 3             F_elk  7667.352
#> 4            Fc_elk  1017.588
#> 5          Findings  5302.763
#> 6 Findings_reindeer  6133.989
#> 7      Findings_elk  8959.729

(sigmaval <- median(opt_sigma_vals$opt_sig))
#> [1] 6133.989
resolu <- 250

Create density raster

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

sh <- st_read("../data/derived_data/vector/sh_nomarsh.shp")
#> Reading layer `sh_nomarsh' from data source `E:\Scales of Transformation\GitLab\b1\data\derived_data\vector\sh_nomarsh.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 2 features and 1 field
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: 489122.6 ymin: 5913462 xmax: 650128.7 ymax: 6084419
#> epsg (SRID):    NA
#> proj4string:    +proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs
cropstack <- mask(dens_stack, sh)

Rescale data by range

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

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