Raster

library(raster)
library(sf)
library(ggplot2)
library(tmap)
library(scales)
library(gridExtra)

Load data

Rasters

Load prepared rasters.

load("1_data/manual/slope_clip.rData")
load("1_data/manual/aspect_clip.rData")
load("1_data/manual/cell_clip.rData")

# slope_clip<-raster("1_data/manual/slope_clip.tif")
# aspect_clip<-raster("1_data/manual/aspect_clip.tif")

Pinch-points

Load buffered pinch-point polygons

# buffered pinch point polygons
load("3_pipeline/store/pp_s_clip_clus_buf_ext.rData")

# study area
load("1_data/manual/study_area.rdata")

# load clipped road data
load("1_data/manual/v_road_clip.rdata")

# load polygon representing North Central Region of Saskatchewan region of the ribbon of green
load("1_data/manual/nscr.rdata")  

# pinch point clusters
load("3_pipeline/tmp/pp_s_clip_clus.rData")

# buffered pinchpoints
load("3_pipeline/store/pp_s_clip_clus_buf.rData")

Prepare rasters

Stack rasters into a single object

#create a list of raster objects to include
r_list<-list(slope_clip, aspect_clip)

cell_clip<-stack(r_list)

#save stack as r object
stackSave(r_stack,file="3_pipeline/store/r_stack")

Extract data

Check the extreme values if the raster stack to identity any values that shouldn’t be included (e.g. -9999).

stackOpen("3_pipeline/store/r_stack")

minValue(r_stack)
maxValue(r_stack)

hist(r_stack)

Summary table

Extract and summarize raster data within each pinch-point polygon

#calculate mean
pp_r_mean<-extract(r_stack,  pp_s_clip_clus_buf_ext, fun=mean, na.rm=TRUE)
colnames(pp_r_mean)<-paste("mean", substr(colnames(pp_r_mean), 1, nchar(colnames(pp_r_mean))-5), sep="_")

#calculate sd
pp_r_sd<-extract(r_stack,  pp_s_clip_clus_buf_ext, fun=sd, na.rm=TRUE)                   
colnames(pp_r_sd)<-paste("sd", substr(colnames(pp_r_sd), 1, nchar(colnames(pp_r_sd))-5), sep="_")          

#calculate max
pp_r_max<-extract(r_stack,  pp_s_clip_clus_buf_ext, fun=max, na.rm=TRUE)                   
colnames(pp_r_max)<-paste("max", substr(colnames(pp_r_max), 1, nchar(colnames(pp_r_sd))-3), sep="_") 

#calculate min
pp_r_min<-extract(r_stack,  pp_s_clip_clus_buf_ext, fun=min, na.rm=TRUE)                   
colnames(pp_r_min)<-paste("min", substr(colnames(pp_r_min), 1, nchar(colnames(pp_r_sd))-3), sep="_")


##### Cell data #########
pp_r_l_mean<-extract(cell_clip,  pp_s_clip_clus_buf_ext, fun=mean, na.rm=TRUE)
colnames(pp_r_l_mean)<-paste("mean_cell", substr(colnames(pp_r_l_mean), 1, nchar(colnames(pp_r_l_mean))-5))

#calculate sd
pp_r_l_sd<-extract(cell_clip,  pp_s_clip_clus_buf_ext, fun=sd, na.rm=TRUE)                   
colnames(pp_r_l_sd)<-paste("sd_cell", substr(colnames(pp_r_l_sd), 1, nchar(colnames(pp_r_l_sd))-5))          

#calculate max
pp_r_l_max<-extract(cell_clip,  pp_s_clip_clus_buf_ext, fun=max, na.rm=TRUE)                   
colnames(pp_r_l_max)<-paste("max_cell", substr(colnames(pp_r_l_max), 1, nchar(colnames(pp_r_l_sd))-3)) 

#calculate min
pp_r_l_min<-extract(cell_clip,  pp_s_clip_clus_buf_ext, fun=min, na.rm=TRUE)                   
colnames(pp_r_l_min)<-paste("min_cell", substr(colnames(pp_r_l_min), 1, nchar(colnames(pp_r_l_sd))-3))



#combine metrics into a single table
pp_r<-cbind(ID=pp_s_clip_clus_buf_ext$ID, pp_r_min, pp_r_max, pp_r_mean, pp_r_sd, pp_r_l_min, pp_r_l_max, pp_r_l_mean, pp_r_l_sd)


pp_r<-cbind(pp_r, pp_r_l_min, pp_r_l_max, pp_r_l_mean, pp_r_l_sd )
#reorder columns
pp_r<-pp_r[,c("ID","min_slope", "max_slope", "mean_slope", "sd_slope", "min_aspect", "max_aspect", "mean_aspect", "sd_aspect", "min_cell", "max_cell", "mean_cell", "sd_cell")]

# save
save(pp_r,file="3_pipeline/store/pp_r.rData")
ID min_slope min_aspect max_slope max_aspect mean_slope mean_aspect sd_slope sd_aspect min_cell max_cell mean_cell sd_cell
1 0.0812411 9.523023 29.51549 350.7695 9.270331 259.04827 7.952051 57.05700 0.0000456 0.0056719 0.0019634 0.0015081
2 0.0704860 3.451188 36.16954 350.3948 9.332548 116.92821 8.964192 117.51785 0.0000672 0.0058609 0.0015163 0.0012145
3 0.0015890 -1.000000 36.34338 355.8282 10.351893 61.10899 8.540351 74.69963 0.0006340 0.0118666 0.0041836 0.0021918
4 0.0302920 7.249906 32.39600 355.5254 10.835116 274.18811 9.079490 99.14873 0.0001906 0.0114184 0.0038657 0.0031939
5 0.5257717 17.686188 15.55570 304.5115 6.774683 71.43680 3.858803 61.59105 0.0002890 0.0057767 0.0022719 0.0014625
6 0.3955552 37.618832 14.68072 347.2162 5.808701 137.50692 3.444783 61.99655 0.0003289 0.0172053 0.0062886 0.0043468
7 2.5336754 158.527191 32.88671 341.4302 14.066558 310.33579 9.399289 26.40259 0.0003269 0.0035523 0.0017525 0.0010484
8 0.1472756 29.808361 29.34639 302.2689 12.175097 119.74625 7.320759 33.05506 0.0003237 0.0137275 0.0033210 0.0027728
9 0.0011139 -1.000000 19.93499 357.6785 5.804861 266.47090 4.539194 88.67824 0.0007211 0.0043416 0.0018871 0.0006696
10 0.5570856 36.425022 11.98081 311.1764 4.675243 82.88598 3.742990 32.59775 0.0001137 0.0002259 0.0001730 0.0000241
11 0.4183576 23.960176 36.29103 316.6076 15.518173 171.71226 7.504589 48.60800 0.0001316 0.0389883 0.0038852 0.0081698
12 0.2824213 8.322483 36.79930 346.5271 12.061293 279.54274 9.658006 34.58862 0.0002528 0.0044315 0.0016714 0.0010472
13 1.4415734 84.780663 29.49452 216.7601 11.417863 176.32245 5.478304 17.11099 0.0009026 0.0554351 0.0064237 0.0077582
14 3.7535818 130.068848 19.73516 173.1658 11.018893 151.23519 3.648061 12.60279 0.0016754 0.0879066 0.0317502 0.0264956

Summary histograms

Create a histogram of raster values for each pinch point.

poly<-pp_s_clip_clus_buf_ext
#poly$ID <- 1:length(poly)
poly$layer <- NULL
d <- data.frame(poly)
v<-extract(r_stack,  poly, na.rm=TRUE, df=TRUE)

vd <- merge(d, v, by="ID")
vd<- vd[-c(3)]
save(vd, file="3_pipeline/tmp/vd.rData")

# x<-vd[vd$ID==1,]
# h<-hist(aspect_clip)

############# Slope plot #################

# basemap
m4<-
  tm_shape(slope_clip)+
    tm_raster(palette = c("white", "#00441b"),
      title = "Slope",
              contrast = c(0,.8),
              style = "cont")+
   tm_shape(v_road_clip)+tm_lines(col="white")+
  tm_layout(bg.col="#d9d9d9")+
  tm_shape(nscr)+
    tm_fill(col="white", alpha=.4)+
    tm_borders(col = "#636363")+
  tm_shape(pp_s_clip_clus)  +
    tm_fill(col="MAP_COLORS", palette = "Dark2", alpha=.8)+
    tm_borders(alpha=1)+
   tm_shape(pp_s_clip_clus_buf)+
    tm_text(text = "group",
          col="black",
          size = 1.3)+
  tm_layout(legend.outside = FALSE, 
            legend.text.size = .8, 
            legend.bg.color="white", 
            legend.frame=TRUE,
            legend.title.size=1, 
            legend.frame.lwd=.8,
            )
m4

#convert tmap to grob
g1<-tmap_grob(m4)


# create hist plots in ggplot convert to grob objects which will be plotted in tmap.

g2<-vd%>%
    ggplot(aes(x = slope_clip)) +
    geom_histogram(binwidth=1, position = "identity", aes(y=..ncount..), fill="#00441b") +
  #geom_density(alpha=.2, fill="#FF6666")+
    scale_y_continuous(labels = percent_format())+
  labs(x="slope", y="% count")+
    theme(#panel.background = element_blank(),
          axis.ticks=element_blank(),
          strip.background =element_rect(fill="white"),
          panel.background = element_rect(fill="#ece2f0"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank())+
  facet_grid(~ID)
  # geom_vline(aes(xintercept=grp.mean),
  #           color="blue", linetype="dashed", size=1)+
g2


lay2 <- rbind(c(1),
              c(1),
              c(1),
              c(1),
              c(2))   

m2<-arrangeGrob(g1, g2, layout_matrix =lay2)

ggsave("4_output/maps/summary/raster_slope_hist_map.png", m2, width=15, height=12
       )

############# Aspect plot #################

m5<-
  tm_shape(aspect_clip)+
    tm_raster(palette= c("white", "#023858"),
              title = "Aspect",
              contrast = c(0,.8))+
  tm_shape(v_road_clip)+tm_lines(col="white")+
  tm_layout(bg.col="#d9d9d9")+
  tm_shape(nscr)+
    tm_fill(col="white", alpha=.4)+
    tm_borders(col = "#636363")+
  tm_shape(pp_s_clip_clus)  +
    tm_fill(col="MAP_COLORS", palette = "Dark2", alpha=.8)+
    tm_borders(alpha=1)+
   tm_shape(pp_s_clip_clus_buf)+
    tm_text(text = "group",
          col="black",
          size = 1.3)+
  
  
  tm_layout(legend.outside = FALSE, 
            legend.text.size = .8, 
            legend.bg.color="white", 
            legend.frame=TRUE,
            legend.title.size=1, 
            legend.frame.lwd=.8,
            )
m5

g3<-tmap_grob(m5)


g4<-vd%>%
    ggplot(aes(x = aspect_clip)) +
    geom_histogram(binwidth=1, position = "identity", aes(y=..ncount..), fill="#023858") +
  #geom_density(alpha=.2, fill="#FF6666")+
    scale_y_continuous(labels = percent_format())+
  labs(x="aspect", y="% count")+
    theme(#panel.background = element_blank(),
          axis.ticks=element_blank(),
          strip.background =element_rect(fill="white"),
          panel.background = element_rect(fill="#ece2f0"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank())+
  facet_grid(~ID)
  # geom_vline(aes(xintercept=grp.mean),
  #           color="blue", linetype="dashed", size=1)+
g4


m3<-arrangeGrob(g3, g4, layout_matrix =lay2)

ggsave("4_output/maps/summary/raster_aspect_hist_map.png", m3, width=15, height=12)



######## Cell phone use ##############

m6<-
  #tm_shape(v_road_clip)+tm_lines(col="#636363")+
  #tm_layout(bg.col="#d9d9d9")+
  tm_shape(nscr)+
    tm_fill(col="white", alpha=.4)+
    tm_borders(col = "#636363")+
  tm_shape(cell_clip)+
    tm_raster(palette = c("white", "#f03b20"),
      title = "visitation density",
              contrast = c(0, 1),
              style = "cont",
      alpha = .8)+
  tm_shape(pp_s_clip_clus)  +
    tm_fill(col="MAP_COLORS", palette = "Dark2", alpha=.8)+
    tm_borders(alpha=1)+
   tm_shape(pp_s_clip_clus_buf)+
    tm_text(text = "group",
          col="black",
          size = 1.3)+
  tm_layout(legend.outside = FALSE, 
            legend.text.size = .8, 
            legend.bg.color="white", 
            legend.frame=TRUE,
            legend.title.size=1, 
            legend.frame.lwd=1,
            )
m6

g5<-tmap_grob(m6)

poly<-pp_s_clip_clus_buf_ext
#poly$ID <- 1:length(poly)
poly$layer <- NULL
d <- data.frame(poly)
v<-extract(cell_clip,  poly, na.rm=TRUE, df=TRUE)

vd <- merge(d, v, by="ID")
vd<- vd[-c(3)]


g6<-vd%>%
    ggplot(aes(x = Visitation_Density_InclRoads_20210625)) +
    geom_histogram(binwidth=.001, position = "identity", aes(y=..ncount..), fill="#023858") +
  #geom_density(alpha=.2, fill="#FF6666")+
    scale_y_continuous(labels = percent_format())+
  labs(x="density of use", y="% count")+
    theme(#panel.background = element_blank(),
          axis.ticks=element_blank(),
          strip.background =element_rect(fill="white"),
          panel.background = element_rect(fill="#ece2f0"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank())+
  facet_grid(~ID)
  # geom_vline(aes(xintercept=grp.mean),
  #           color="blue", linetype="dashed", size=1)+
g6


m7<-arrangeGrob(g5, g6, layout_matrix =lay2)

ggsave("4_output/maps/summary/raster_cell_hist_map.png", m7, width=15, height=12)
Histograms of slope values for each pinch point.

Figure 3.8: Histograms of slope values for each pinch point.

Histograms of aspect values for each pinch point.

Figure 3.9: Histograms of aspect values for each pinch point.

Histograms of cellphone density values for each pinch point.

Figure 3.10: Histograms of cellphone density values for each pinch point.