library(magrittr)
library(sf)
library(raster)
library(FuzzyLandscapes)
library(dplyr)
library(ggplot2)
library(spatstat)
library(maptools)
library(mapview)

How many points should be selected in which radius

n_points <- 100
radius <- 2000
n <- 4 

Reindeer Hunters

Random Sampling

fr_wet_depression_near_flint <- raster("../data/derived_data/raster/FuzzyPot/wet_depression_near_flint.tif")

set.seed(1234)
fitloc <- Map(
  function(x){
    mod_fr_reindeer <- fr_wet_depression_near_flint/(18794919/(1+x))
    mod_fr_reindeer[is.na(mod_fr_reindeer[])] <- 0 
    
    set.seed(1234)
    rpoints_reindeer <- rpoispp(lambda = maptools::as.im.RasterLayer(mod_fr_reindeer))
    
    frp_reindeer <- as.SpatialPoints.ppp(rpoints_reindeer) 
    proj4string(frp_reindeer) <- proj4string(fr_wet_depression_near_flint)
    return(data.frame(x=c(x),
                      np = length(frp_reindeer)))
  },
  x=c(1:2))


fitloc <- do.call("rbind",fitloc)
x <- predict(lm(x~np,fitloc),data.frame(np=n_points))


mod_fr_reindeer <- fr_wet_depression_near_flint/(18794919/(1+x))
mod_fr_reindeer[is.na(mod_fr_reindeer[])] <- 0 
set.seed(1234)
rpoints_reindeer <- rpoispp(lambda = maptools::as.im.RasterLayer(mod_fr_reindeer))

fp_reindeer <- as.SpatialPoints.ppp(rpoints_reindeer)
proj4string(fp_reindeer) <- proj4string(fr_wet_depression_near_flint)
length(fp_reindeer)
#> [1] 100

plot(fr_wet_depression_near_flint);points(fp_reindeer,pch=16,cex=0.5)


F_reindeer <- rpoints_reindeer
F_reindeer$marks <- as.factor(rep("Random_reindeer",F_reindeer$n))

Random Clustered Poisson Sampling

fitloc <- Map(
  function(x){
    mod_fr_reindeer <- fr_wet_depression_near_flint/(18794919/(1+x))
    mod_fr_reindeer[is.na(mod_fr_reindeer[])] <- 0 
    
    set.seed(1234)
    rpc_reindeer <- rPoissonCluster(kappa = maptools::as.im.RasterLayer(mod_fr_reindeer),
                                    expand = 0.1,
                                    win = trim.rectangle(as.rectangle(as.owin(maptools::as.im.RasterLayer(mod_fr_reindeer))),100),
                                    rcluster = function(x0, y0, radius, n) {
                                      return(runifdisc(n, radius, centre=c(x0, y0)))
                                    }, 
                                    radius=radius, 
                                    n=n)
    
    frp_reindeer <- as.SpatialPoints.ppp(rpc_reindeer)
    proj4string(frp_reindeer) <- proj4string(fr_wet_depression_near_flint)
    return(data.frame(x=c(x),
                      np = length(frp_reindeer)))
  },
  x=c(1:2))


fitloc <- do.call("rbind",fitloc)
x <- predict(lm(x~np,fitloc),data.frame(np=n_points))


mod_fr_reindeer <- fr_wet_depression_near_flint/(18794919/(1+x))
mod_fr_reindeer[is.na(mod_fr_reindeer[])] <- 0 

set.seed(1234)
rpc_reindeer <- rPoissonCluster(kappa = maptools::as.im.RasterLayer(mod_fr_reindeer),
                                expand = 0.1,
                                win = trim.rectangle(as.rectangle(as.owin(maptools::as.im.RasterLayer(mod_fr_reindeer))),100),
                                rcluster = function(x0, y0, radius, n) {
                                  return(runifdisc(n, radius, centre=c(x0, y0)))
                                }, 
                                radius=radius, 
                                n=n)

frp_reindeer <- as.SpatialPoints.ppp(rpc_reindeer)
proj4string(frp_reindeer) <- proj4string(fr_wet_depression_near_flint)
length(frp_reindeer)
#> [1] 112

plot(fr_wet_depression_near_flint);points(frp_reindeer,pch=16,cex=0.5)


Fc_reindeer <- rpc_reindeer
Fc_reindeer$marks <- as.factor(rep("Cluster_reindeer",Fc_reindeer$n))
saveRDS(Fc_reindeer, file = "../data/derived_data/vector/Fc_reindeer.rds")

F_reindeer$window <- Fc_reindeer$window
saveRDS(F_reindeer, file = "../data/derived_data/vector/F_reindeer.rds")

mapview(fp_reindeer)+mapview(frp_reindeer)

Elk Hunters

Random Sampling

fr_wet_sandy_ridge_near_flint <- raster("../data/derived_data/raster/FuzzyPot/wet_sandy_ridge_near_flint.tif")

set.seed(1234)
fitloc <- Map(
  function(x){
    mod_fr_elk <- fr_wet_sandy_ridge_near_flint/(18794919/(1+x))
    mod_fr_elk[is.na(mod_fr_elk[])] <- 0 
    
    set.seed(1234)
    rpoints_elk <- rpoispp(lambda = maptools::as.im.RasterLayer(mod_fr_elk))
    
    frp_elk <- as.SpatialPoints.ppp(rpoints_elk) 
    proj4string(frp_elk) <- proj4string(fr_wet_sandy_ridge_near_flint)
    return(data.frame(x=c(x),
                      np = length(frp_elk)))
  },
  x=c(1:2))


fitloc <- do.call("rbind",fitloc)
x <- predict(lm(x~np,fitloc),data.frame(np=n_points))


mod_fr_elk <- fr_wet_sandy_ridge_near_flint/(18794919/(1+x))
mod_fr_elk[is.na(mod_fr_elk[])] <- 0
set.seed(1234)
rpoints_elk <- rpoispp(lambda = maptools::as.im.RasterLayer(mod_fr_elk))


fp_elk <- as.SpatialPoints.ppp(rpoints_elk)
proj4string(fp_elk) <- proj4string(fr_wet_sandy_ridge_near_flint)
length(fp_elk)
#> [1] 100

plot(fr_wet_sandy_ridge_near_flint);points(fp_elk,pch=16,cex=0.5)


F_elk <- rpoints_elk
F_elk$marks <- as.factor(rep("Random_elk",F_elk$n))

F_elk$window <- Fc_reindeer$window
saveRDS(F_elk, file = "../data/derived_data/vector/F_elk.rds")

Random Clustered Poisson Sampling

fitloc <- Map(
  function(x){
    mod_fr_elk <- fr_wet_sandy_ridge_near_flint/(18794919/(1+x))
    mod_fr_elk[is.na(mod_fr_elk[])] <- 0 
    
    set.seed(1234)
    rpc_elk <- rPoissonCluster(kappa = maptools::as.im.RasterLayer(mod_fr_elk),
                               expand = 0.1,
                               win = trim.rectangle(as.rectangle(as.owin(maptools::as.im.RasterLayer(mod_fr_elk))),100),
                               rcluster = function(x0, y0, radius, n) {
                                 return(runifdisc(n, radius, centre=c(x0, y0)))
                               }, 
                               radius=radius, 
                               n=n)
    
    frp_elk <- as.SpatialPoints.ppp(rpc_elk)
    proj4string(frp_elk) <- proj4string(fr_wet_sandy_ridge_near_flint)
    return(data.frame(x=c(x),
                      np = length(frp_elk)))
  },
  x=c(1:2))


fitloc <- do.call("rbind",fitloc)
x <- predict(lm(x~np,fitloc),data.frame(np=n_points))


mod_fr_elk <- fr_wet_sandy_ridge_near_flint/(18794919/(1+x))
mod_fr_elk[is.na(mod_fr_elk[])] <- 0 

set.seed(1234)
rpc_elk <- rPoissonCluster(kappa = maptools::as.im.RasterLayer(mod_fr_elk),
                           expand = 0.1,
                           win = trim.rectangle(as.rectangle(as.owin(maptools::as.im.RasterLayer(mod_fr_elk))),100),
                           rcluster = function(x0, y0, radius, n) {
                             return(runifdisc(n, radius, centre=c(x0, y0)))
                           }, 
                           radius=radius, 
                           n=n)

frp_elk <- as.SpatialPoints.ppp(rpc_elk)
proj4string(frp_elk) <- proj4string(fr_wet_sandy_ridge_near_flint)
length(frp_elk)
#> [1] 100

plot(fr_wet_sandy_ridge_near_flint);points(frp_elk,pch=16,cex=0.5)


Fc_elk <- rpc_elk
Fc_elk$marks <- as.factor(rep("Cluster_elk",Fc_elk$n))
saveRDS(Fc_elk, file = "../data/derived_data/vector/Fc_elk.rds")

mapview(fp_elk)+mapview(frp_elk)

Findings

b1_locations_centroid <- st_read("../data/derived_data/vector/b1_locations_centroid.gpkg", 
                                 "b1_locations_centroid")%>% 
  .[-which(is.na(raster::extract(fr_wet_sandy_ridge_near_flint,.))),]
#> Reading layer `b1_locations_centroid' from data source `E:\Scales of Transformation\GitLab\b1\data\derived_data\vector\b1_locations_centroid.gpkg' using driver `GPKG'
#> Simple feature collection with 68 features and 1 field
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 8.370252 ymin: 53.45151 xmax: 10.92449 ymax: 54.88782
#> epsg (SRID):    4326
#> proj4string:    +proj=longlat +datum=WGS84 +no_defs
ppp.find <- function(dat){
  Findings <- dplyr::select(dat,1) %>% 
    as(.,'Spatial') %>% 
    spTransform(., CRSobj = CRS(proj4string(fr_wet_depression_near_flint)))
  
  Findings$la_id <- factor("Findings")
  Findings <- as(Findings, "ppp")
  Findings$window <- Fc_reindeer$window
  return(Findings)
}


Findings <- ppp.find(b1_locations_centroid)
saveRDS(Findings, file = "../data/derived_data/vector/Findings.rds")

Findings_reindeer <- ppp.find(b1_locations_centroid %>% dplyr::filter(Time == "Reindeerhunter"))
saveRDS(Findings_reindeer, file = "../data/derived_data/vector/Findings_reindeer.rds")

fire_kest <- Kest(Findings_reindeer)
fire_fest <- Fest(Findings_reindeer)
fire_gest <- Gest(Findings_reindeer)

Findings_elk <- ppp.find(b1_locations_centroid %>% dplyr::filter(Time == "Elkhunter"))
saveRDS(Findings_elk, file = "../data/derived_data/vector/Findings_elk.rds")

fielk_kest <- Kest(Findings_elk)
fielk_fest <- Fest(Findings_elk)
fielk_gest <- Gest(Findings_elk)

ggplot(data=rbind(Kest(Findings) %>% mutate(Type = "Both"),
                  fire_kest %>% mutate(Type = "Reindeeer"),
                  fielk_kest %>% mutate(Type = "Elkhunters"),
                  fire_kest %>% mutate(Type = "Theoretical",
                                       border = theo,
                                       trans = theo,
                                       iso = theo)),
       aes(x=r, y=trans, colour=Type)) + 
  geom_line() +
  theme_bw() +
  coord_cartesian(xlim = c(0,20000),ylim=c(0,4e+09))+
  labs(y = "K(r)")

Compare based upon Kest Function

compda <- do.call("rbind",Map(function(x,ty,va,co){Kest(get(x)) %>% mutate(Type = ty,
                                                                           Variant = va,
                                                                           Comp = co)},
                              x = c("Findings_reindeer","Findings_reindeer",
                                    "Findings_elk","Findings_elk",
                                    "F_reindeer","Fc_reindeer",
                                    "F_elk","Fc_elk"),
                              ty = c("Findings","Findings",
                                     "Findings","Findings",
                                     "Fuzzy Reindeer hunters","Fuzzy Reindeer hunters",
                                     "Fuzzy Elk hunters","Fuzzy Elk hunters"),
                              va = c("Non-clustered","Clustered",
                                     "Non-clustered","Clustered",
                                     "Non-clustered","Clustered",
                                     "Non-clustered","Clustered"),
                              co = c("Reindeerhunters","Reindeerhunters",
                                     "Elkhunters","Elkhunters",
                                     "Reindeerhunters","Reindeerhunters",
                                     "Elkhunters","Elkhunters"
                              )
))

compda_r <- rbind(compda,
                  Kest(Findings_reindeer) %>% 
                    mutate(Type = "Theoretical",
                           border = theo,
                           trans = theo,
                           iso = theo,
                           Variant = "Non-clustered",
                           Comp = "Reindeerhunters"),
                  Kest(Findings_reindeer) %>% 
                    mutate(Type = "Theoretical",
                           border = theo,
                           trans = theo,
                           iso = theo,
                           Variant = "Clustered",
                           Comp = "Reindeerhunters"),
                  Kest(Findings_reindeer) %>% 
                    mutate(Type = "Theoretical",
                           border = theo,
                           trans = theo,
                           iso = theo,
                           Variant = "Non-clustered",
                           Comp = "Elkhunters"),
                  Kest(Findings_reindeer) %>% 
                    mutate(Type = "Theoretical",
                           border = theo,
                           trans = theo,
                           iso = theo,
                           Variant = "Clustered",
                           Comp = "Elkhunters"))

plot_ripley <- ggplot(data=compda_r %>% filter(Type != "Binary addition"), aes(x=r, y=trans, colour=Type)) + 
  geom_line(aes(linetype=Type, color=Type)) +
  scale_linetype_manual(values=c("dashed","solid","solid","dotted")) +
  theme_bw() +
  coord_cartesian(xlim = c(0,40000),ylim=c(0,10e+09))+
  labs(y = "K(r)")+
  facet_grid(rows = vars(Comp,Variant))
plot_ripley

Compare based upon Gest Function

compda <- do.call("rbind",Map(function(x,ty,va,co){Gest(get(x)) %>% mutate(Type = ty,
                                                                           Variant = va,
                                                                           Comp = co)},
                              x = c("Findings_reindeer","Findings_reindeer",
                                    "Findings_elk","Findings_elk",
                                    "F_reindeer","Fc_reindeer",
                                    "F_elk","Fc_elk"),
                              ty = c("Findings","Findings",
                                     "Findings","Findings",
                                     "Fuzzy Reindeer hunters","Fuzzy Reindeer hunters",
                                     "Fuzzy Elk hunters","Fuzzy Elk hunters"),
                              va = c("Non-clustered","Clustered",
                                     "Non-clustered","Clustered",
                                     "Non-clustered","Clustered",
                                     "Non-clustered","Clustered"),
                              co = c("Reindeerhunters","Reindeerhunters",
                                     "Elkhunters","Elkhunters",
                                     "Reindeerhunters","Reindeerhunters",
                                     "Elkhunters","Elkhunters"
                              )
))


compda_g <- rbind(compda,
                  Gest(Findings_reindeer) %>% 
                    mutate(Type = "Theoretical",
                           han = theo,
                           rs = theo,
                           km = theo,
                           Variant = "Non-clustered",
                           Comp = "Reindeerhunters"),
                  Gest(Findings_reindeer) %>% 
                    mutate(Type = "Theoretical",
                           han = theo,
                           rs = theo,
                           km = theo,
                           Variant = "Clustered",
                           Comp = "Reindeerhunters"),
                  Gest(Findings_reindeer) %>% 
                    mutate(Type = "Theoretical",
                           han = theo,
                           rs = theo,
                           km = theo,
                           Variant = "Non-clustered",
                           Comp = "Elkhunters"),
                  Gest(Findings_reindeer) %>% 
                    mutate(Type = "Theoretical",
                           han = theo,
                           rs = theo,
                           km = theo,
                           Variant = "Clustered",
                           Comp = "Elkhunters"))

plot_g <- ggplot(data=compda_g %>% filter(Type != "Binary addition"), aes(x=r, y=km, colour=Type)) + 
  geom_line(aes(linetype=Type, color=Type)) +
  scale_linetype_manual(values=c("dashed","solid","solid","dotted")) +
  theme_bw() +
  coord_cartesian(xlim = c(0,20000),ylim=c(0,1))+
  labs(y = "G(r)")+
  facet_grid(rows = vars(Comp,Variant))
plot_g

Compare based upon Fest Function

compda <- do.call("rbind",Map(function(x,ty,va,co){Fest(get(x)) %>% mutate(Type = ty,
                                                                           Variant = va,
                                                                           Comp = co)},
                              x = c("Findings_reindeer","Findings_reindeer",
                                    "Findings_elk","Findings_elk",
                                    "F_reindeer","Fc_reindeer",
                                    "F_elk","Fc_elk"),
                              ty = c("Findings","Findings",
                                     "Findings","Findings",
                                     "Fuzzy Reindeer hunters","Fuzzy Reindeer hunters",
                                     "Fuzzy Elk hunters","Fuzzy Elk hunters"),
                              va = c("Non-clustered","Clustered",
                                     "Non-clustered","Clustered",
                                     "Non-clustered","Clustered",
                                     "Non-clustered","Clustered"),
                              co = c("Reindeerhunters","Reindeerhunters",
                                     "Elkhunters","Elkhunters",
                                     "Reindeerhunters","Reindeerhunters",
                                     "Elkhunters","Elkhunters"
                              )
))


compda_f <- rbind(compda,
                  Fest(Findings_reindeer) %>% 
                    mutate(Type = "Theoretical",
                           cs = theo,
                           rs = theo,
                           km = theo,
                           Variant = "Non-clustered",
                           Comp = "Reindeerhunters"),
                  Fest(Findings_reindeer) %>% 
                    mutate(Type = "Theoretical",
                           cs = theo,
                           rs = theo,
                           km = theo,
                           Variant = "Clustered",
                           Comp = "Reindeerhunters"),
                  Fest(Findings_reindeer) %>% 
                    mutate(Type = "Theoretical",
                           cs = theo,
                           rs = theo,
                           km = theo,
                           Variant = "Non-clustered",
                           Comp = "Elkhunters"),
                  Fest(Findings_reindeer) %>% 
                    mutate(Type = "Theoretical",
                           cs = theo,
                           rs = theo,
                           km = theo,
                           Variant = "Clustered",
                           Comp = "Elkhunters"))

plot_f <- ggplot(data=compda_f %>% filter(Type != "Binary addition"), aes(x=r, y=cs, colour=Type)) + 
  geom_line(aes(linetype=Type, color=Type)) +
  scale_linetype_manual(values=c("dashed","solid","solid","dotted")) +
  theme_bw() +
  coord_cartesian(xlim = c(0,40000),ylim=c(0,1))+
  labs(y = "F(r)")+
  facet_grid(rows = vars(Comp,Variant))
plot_f

Combined plots

library(ggpubr)
ggsave("../plots/plot_kgf.png", 
       plot = ggarrange(plot_g+theme(legend.title=element_blank())+ 
                          theme(
                            strip.background = element_blank(),
                            strip.text.y = element_blank()
                          ),
                        plot_f + theme(legend.position = "none")+ 
                          theme(
                            strip.background = element_blank(),
                            strip.text.y = element_blank()
                          ), 
                        plot_ripley + theme(legend.position = "none"),
                        labels = c("G", "F", "K"),
                        ncol = 3, nrow = 1,
                        legend="bottom",
                        common.legend = TRUE,
                        widths = c(1,1, 1.2)),
       width = 9, height = 9)