# Script Information ------------------------------------------------------
##
## Description: code to run analyses and figures for manuscript "Advancements in monitoring: a comparison of traditional and application-based tools for measuring outdoor recreation"
##
## Author: Talia Vilalta Capdevila, Research Technician, Yellowstone to Yukon Conservation Initiative, MSc 
## 
## Date last modified: February 20, 2024
##
## Manuscript corresponding author: Brynn McLellan, brynn@y2y.net, Research Associate, Yellowstone to Yukon Conservation Initiative, MSc  
## 
## Notes: 
##
## The trail camera data, trail counter data, and Strava Metro data are not openly available due to sensitivity and license agreements from data contributors. 
## Data are, however, available from data contributors upon request and possible data license agreements. 
## Please see data availability section of this publication for the contact information of data contributors. 
##  
### description of data set names in this script:
#### camera =  trail camera 
#### TRAFx = infrared trail counter
#### STRAVA or Strava = Strava Metro 
#### Strava Heatmap or heatmap = Strava Global Heatmap 


# Setting up the Environment --------------------------------------------------

# clear global environment
rm(list=ls())

# loading packages 
library(broom) # turn into tidy tibbles
library(dplyr) # for using tidyverse
library(data.table) # for data tables 
library(ggplot2) # for graphs 
library(gridExtra) # for graphs 
library(hrbrthemes) # additional themes and components for ggplot 
library(lubridate) # manipulating dates
library(patchwork) # for graphs
library(plotly) # for interactive maps
library(plyr) # to count levels
library(scales) # for labeling x axes with commas
library(tidyr) # for data wrangling
library(tidyverse) # to use tidyverse
library(ggpubr) # for plotting 
library(zoo) # S3 class with methods for indexed totally ordered observations

# set working directory to folder where RScript is saved
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))


# Load in Required Data ---------------------------------------------------

# load in spatial matches of tools 
## layers were generated in ArcMap and represent locations where tools were within a distance of each other to compare counts from pairs of tools 
## 30m buffer to match Strava Metro trails to cameras/counters, 200m to match cameras and counters

# camera : TRAFx
Cam_TRAFx_200m_Matches <- read.csv("./Spatial_Matches/Cam_TRAFx_200m_Matches.csv")

# camera : STRAVA
## for STRAVA: some locations have more than 1 edgeuid due to Strava Metro map updates over the three year of data (2017, 2018, 2019)
## edgeuid is the unique identification number each Strava Metro segment has, as Strava Metro updates with new OpenStreetMap data, all the edgeuids change
Cam_STRAVA_Matches_30m <- read.csv("./Spatial_Matches/Camera_STRAVA_Matches_30m.csv")

# TRAFx : STRAVA
TRAFx_STRAVA_Matches_30m <- read.csv("./Spatial_Matches/TRAFx_STRAVA_Matches_30m.csv")


# Creating Unique Data Sets for Each Activity and Tool ----------------

#### TRAFx data wrangling ####

# TRAFx daily count data 
## load TRAFx daily count data; select location, date and count columns, change date format and put location names to lowercase
TRAFx_Data <- read.csv("./Count_Data_MS/All_TRAFx_Counts_Daily_20230525.csv") %>%
dplyr::select(Location,Daily_Count, date)%>%
  mutate(date = as.POSIXct(date, format = "%Y-%m-%d"))%>%
  mutate(Location = tolower(Location))
## total number of TRAFx counts
sum(TRAFx_Data$Daily_Count)

# TRAFx biking counts
## these locations were identified as being bike counters
TRAFX_Biking <-TRAFx_Data %>%
  dplyr::filter(Location %like% 'vc|mb|vehicle')

# TRAFX pedestrian counts 
## these locations were identified as being pedestrian counters (all that were not bike counters)
TRAFx_Pedestrian <- TRAFx_Data %>%
  dplyr::filter(!Location%like% 'vc|mb|vehicle')


#### STRAVA data Wrangling ####

# load STRAVA daily count data 
STRAVA_Data <- read.csv("./Count_Data_MS/Daily_STRAVA_Data20230412.csv") %>%
  mutate(date = as.POSIXct(date, format = "%Y-%m-%d"))
## number of unique Strava segments used
n_distinct(STRAVA_Data$edge_uid)
## number of activities recorded
sum(STRAVA_Data$Daily_Tot_Count)


#### Camera data wrangling ####

# load camera daily count data
Camera_Data_1 <- read.csv("./Count_Data_MS/CameraData_Daily_PCA_ABPARKS_CovInt.csv") 

# camera all activities counts data 
## select columns with relevant information, group by location and date and add up total counts
Camera_Data <- Camera_Data_1 %>%
  dplyr::select(location, dateimage, totaldaily, species) %>% # species refers to the type of human activity detected
  group_by(location, dateimage) %>%
  dplyr::summarise(TotalDaily = sum(c(totaldaily), na.rm = TRUE)) %>%
  ungroup() %>%
  separate(dateimage, c("year", "month","day"),
           remove = FALSE) %>%
  mutate(dateimage = as.POSIXct(dateimage, format = "%Y-%m-%d"))
## total number of recreation activities collected by cameras
sum(Camera_Data$TotalDaily)

# camera biking counts 
## create a data frame of "species" to view all categories
Cam_Species <- as.data.frame(unique(Camera_Data_1$species))
## bike data was considered to be if image had a biker, select only those species 
Camera_Data_Biking <- Camera_Data_1 %>%
  dplyr::filter(species %like% 'biker| e biker| fatbiker')
## select columns, group and add up counts by date and location
Camera_Data_Biking_1 <- Camera_Data_Biking %>%
  dplyr::select(location, dateimage, totaldaily, species) %>%
  dplyr::group_by(location,dateimage) %>%
  dplyr::summarise(TotalDaily = sum(c(totaldaily), na.rm = TRUE)) %>%
  dplyr::ungroup() %>%
  separate(dateimage, c("year", "month","day"),
           remove = FALSE) %>%
  dplyr::mutate(dateimage = as.POSIXct(dateimage, format = "%Y-%m-%d"))


# camera pedestrian counts
## pedestrian data was considered to be if image had a a person on foot, select only those species 
Camera_Data_Pedestrian <- Camera_Data_1 %>%
  dplyr::filter(species %like% 'angler|climber|hiker| hunter|runner|skier|skijorer|snowshoer')
## select columns, group and add up counts by date and location
Camera_Data_Pedestrian_1 <- Camera_Data_Pedestrian %>%
  dplyr::select(location, dateimage, totaldaily, species) %>%
  dplyr::group_by(location, dateimage) %>%
  dplyr::summarise(TotalDaily = sum(c(totaldaily), na.rm = TRUE)) %>%
  dplyr:: ungroup()%>%
  separate(dateimage, c("year", "month","day"),
           remove = FALSE) %>%
  dplyr::mutate(dateimage = as.POSIXct(dateimage, format = "%Y-%m-%d"))



# Joining Data to Spatial Matches -----------------------------------------

#### Camera and TRAFx data joining ####

# all activities 
## joining camera data to spatial match
Cam_200m_Data <- full_join(Cam_TRAFx_200m_Matches, Camera_Data, by = c("Cam_Location" = "location")) %>%
  dplyr::rename(Cam_Daily = TotalDaily, date = dateimage) 

## joining TRAFx data to spatial match
Cam_TRAFx_200m_Data <- full_join(Cam_200m_Data, TRAFx_Data, by = c("Location", "date")) %>%
  dplyr::rename(TRAFx_Daily = Daily_Count, TRAFx_Location = Location) %>%
  dplyr:: mutate(Weekday = weekdays(date)) %>%
  dplyr::filter(Cam_Location != "cpchrzlvelegreii4 x46z42g48 x4666z4g42") %>% #removing match (duplicate site/day)
  dplyr::filter(!is.na(TRAFx_Daily)) %>%
  dplyr::filter(!is.na(Distance))

## combining entries where camera was matched to a VC (vehicle) and TC (pedestrian) counter
## only summing TRAFx values (by camera and date) when they contain vc or ve in  name, otherwise put count values as 99999 - in a new column. Original rows still present and will have to be removed later.
Cam_TRAFx_200m_Data1 <- Cam_TRAFx_200m_Data %>% 
  group_by(Cam_Location, date) %>% 
  dplyr::mutate(TRAFx_Combined = if_else(( TRAFx_Location %like% 'vc|ve')
                                         , sum(TRAFx_Daily)
                                         , as.integer(99999))) %>% ungroup()

## keeping smallest distance when TRAFx are matched to multiple cameras on the same day (there are duplicate camera locations on the same day)
## first arrange by distance - filter will keep the first of the duplicates, then arrange TRAFx location in descending order so those so be will be chosen (has summed counts)
Cam_TRAFx_200m_Data2 <- Cam_TRAFx_200m_Data1 %>%
  dplyr::arrange(date, Distance, desc(TRAFx_Location))

Cam_TRAFx_200m_Data2 <- Cam_TRAFx_200m_Data2 %>%
  dplyr::filter(duplicated(Cam_TRAFx_200m_Data2[, c(3, 11)]) == FALSE)  %>% # remove duplicated camera (column 3)/date (column 11) combos - keeps the first entry
  dplyr::mutate(TRAFx_Fixed = if_else((TRAFx_Combined == '99999')
                                      , TRAFx_Daily
                                      , TRAFx_Combined)) # add a new column that keeps the original counts (tagged as 99999) and puts the combined counts where there were duplicates.


## checking for TRAFx matches to multiple cameras
TRAFx_MultiCam2 <- duplicated(Cam_TRAFx_200m_Data2[, c(5, 11)]) # set columns to match TRAFx location and date
TRAFx_MultiCam2 <- Cam_TRAFx_200m_Data2[TRAFx_MultiCam2, ]

## checking for camera matches to multiple TRAFx
Cam_MultiTRAFx2 <- duplicated(Cam_TRAFx_200m_Data2[,c(3,11)]) # set columns to match camera locations and date
Cam_MultiTRAFx2 <- Cam_TRAFx_200m_Data2[Cam_MultiTRAFx2, ]

## removing 2 cameras over 100m from a double sampled TRAFx (there was a cluster of camera:TRAFx, so initial cleaning still left behind a couple matches)
## removing one entry where a day is duplicated
Cam_TRAFx_200m_Data3 <- Cam_TRAFx_200m_Data2 %>%
  dplyr::filter(!Cam_Location %like% 'bv 22g4 ne-b|bv 22g4 ne-zl') %>%
  dplyr::filter(!Cam_Location %like% 'bv 22g4 ne-b|bv 22g4 ne-zl') %>%
  dplyr::filter(!(Cam_Location =='bv 2x4u8 nw' & year == '2019' & month =='04' & day =='05'))

## number of camera:TRAFx comparisons for all activities
n_distinct(Cam_TRAFx_200m_Data3$Cam_Location) 


# Biking activities 
## join specific camera biking data to spatial matches
Cam_200m_BikingData <- full_join(Cam_TRAFx_200m_Matches, Camera_Data_Biking_1, by = c("Cam_Location" = "location")) %>%
  dplyr::rename(Cam_Daily=TotalDaily, date = dateimage)

## join specific TRAFx biking data, and remove any entries that had NAs in either camera data (used distance column), or TRAFx data (TRAFx_Daily column)
Cam_TRAFx_BikingData <- full_join(Cam_200m_BikingData, TRAFX_Biking, by = c("Location", "date")) %>%
  dplyr::rename(TRAFx_Daily = Daily_Count, TRAFx_Location = Location) %>%
  dplyr::filter(!is.na(TRAFx_Daily)) %>%
  dplyr::filter(!is.na(Distance))

## number of camera:TRAFx comparisons for biking 
n_distinct(Cam_TRAFx_BikingData$Cam_Location)


# Pedestrian activities
## join specific pedestrian camera data to spatial matches
Cam_200m_PedestrianData <- full_join(Cam_TRAFx_200m_Matches, Camera_Data_Pedestrian_1, by = c("Cam_Location" = "location")) %>%
  dplyr::rename(Cam_Daily=TotalDaily, date = dateimage)

## join pedestrian TRAFx data and remove NAs
Cam_TRAFx_PedestrianData <- full_join(Cam_200m_PedestrianData, TRAFx_Pedestrian, by = c("Location",  "date")) %>%
  dplyr::rename(TRAFx_Daily = Daily_Count, TRAFx_Location = Location) %>%
  dplyr::filter(!is.na(TRAFx_Daily)) %>%
  dplyr::filter(!is.na(Distance))

## number of camera TRAFx comparisons for pedestrian
n_distinct(Cam_TRAFx_PedestrianData$Cam_Location)


#### Camera and STRAVA data joining ####

# all activities 
## adding camera data to spatial matches 
Cam_STRAVA_Data1 <- full_join(Cam_STRAVA_Matches_30m, Camera_Data, by = c("Location" = "location")) %>%
  dplyr::rename(Cam_Daily = TotalDaily)

## adding STRAVA data to matches
Cam_STRAVA_Data <- full_join(Cam_STRAVA_Data1, STRAVA_Data, by = c("EdgeUID" = "edge_uid", "dateimage"="date")) %>%
  dplyr::select(-c(Date_Day, Date_Month, Date_Year)) %>%
  dplyr::rename(STRAVA_Daily = Daily_Tot_Count) %>%
  dplyr::filter(!is.na(STRAVA_Daily)) %>% # remove all locations where there was no Strava data 
  filter(!is.na(Cam_Daily))  # remove all locations without camera data 

## remove TCH (Trans Canada highway) matches since cameras were in underpasses, STRAVA is on Legacy trail above - they are not counting the same people
Cam_STRAVA_Data <- Cam_STRAVA_Data %>%
  dplyr::filter(!Location %like% 'tch ') %>%
  mutate(dateimage = as.POSIXct(dateimage, format = "%Y-%m-%d")) %>%
  dplyr::filter(!is.na(STRAVA_Daily)) %>%
  dplyr::filter(!is.na(Cam_Daily))

## number camera : STRAVA comparisons for all activities 
n_distinct(Cam_STRAVA_Data$EdgeUID)


# Biking activities 
## adding unique camera data to spatial matches
Cam_STRAVA_Biking_Data2 <- full_join(Cam_STRAVA_Matches_30m, Camera_Data_Biking_1, by = c("Location" = "location")) %>%
  dplyr::rename(Cam_Bike_Daily = TotalDaily)

## add STRAVA data, and remove any NAs in counts
Cam_STRAVA_Biking_Data3 <- full_join(Cam_STRAVA_Biking_Data2, STRAVA_Data, by = c("EdgeUID" = "edge_uid", "dateimage"="date")) %>%
  dplyr::filter(!is.na(Daily_Ride_Sum)) %>%
  dplyr::filter(!is.na(Cam_Bike_Daily))

## number of camera : STRAVA comparisons for biking 
n_distinct(Cam_STRAVA_Biking_Data3$EdgeUID)


# Pedestrian activities 
## join unique pedestrian Strava data to all camera:Strava data, and remove any entries that do not have Strava pedestrian data
Cam_STRAVA_Pedestrian_Data2 <- full_join(Cam_STRAVA_Matches_30m, Camera_Data_Pedestrian_1, by = c("Location" = "location"))%>%
  dplyr::rename(Cam_Daily = TotalDaily)

## add STRAVA data, and remove any NAs in counts
Cam_STRAVA_Pedestrian_Data3 <- full_join(Cam_STRAVA_Pedestrian_Data2, STRAVA_Data, by = c("EdgeUID" = "edge_uid", "dateimage"="date")) %>%
  dplyr::filter(!is.na(Daily_Ped_Sum)) %>%
  dplyr::filter(!is.na(Cam_Daily))

## number of camera : STRAVA comparisons for pedestrians 
n_distinct(Cam_STRAVA_Pedestrian_Data3$EdgeUID)


#### TRAFx and STRAVA data joining ####

# all activities 
## adding TRAFx data to matches 
TRAFx_STRAVA_Data2 <- full_join(TRAFx_STRAVA_Matches_30m, TRAFx_Data, by = c("Location"))%>%
  dplyr:: rename(TRAFx_Daily = Daily_Count)

## adding Strava data to matches, remove entries with NA values from either tool 
TRAFx_STRAVA_Data3 <- full_join(TRAFx_STRAVA_Data2, STRAVA_Data, by = c("EdgeUID" = "edge_uid", "date"="date")) %>%
  dplyr:: rename(STRAVA_Daily = Daily_Tot_Count) %>%
  filter(!is.na(STRAVA_Daily))%>%
  filter(!is.na(TRAFx_Daily))

## since some vehicle counters and trail counters where in the same Strava area (within 30m of the same segment), their totals have to be combined
TRAFx_STRAVA_Data4 <- TRAFx_STRAVA_Data3 %>% 
  group_by(EdgeUID,date) %>% 
  dplyr::mutate(TRAFx_Combined = if_else((Location %like% 'vc|ve')
                                         , sum(TRAFx_Daily)
                                         , as.integer(99999))) %>% ungroup()

TRAFx_STRAVA_Data5 <- TRAFx_STRAVA_Data4 %>%
  dplyr::arrange(date, desc(Location))

TRAFX_STRAVA_AllData6 <- TRAFx_STRAVA_Data5  %>% 
  dplyr::filter(duplicated(TRAFx_STRAVA_Data5[,c(4,9)]) == FALSE)  %>% # remove if TRAFx location and date are the same
  dplyr::mutate(TRAFx_Fixed = if_else((TRAFx_Combined == '99999')
                                      , TRAFx_Daily
                                      , TRAFx_Combined))

## number of TRAFx : STRAVA comparisons for all activities 
n_distinct(TRAFX_STRAVA_AllData6$EdgeUID)

# biking activities
TRAFx_STRAVA_Biking_Data2 <- full_join(TRAFx_STRAVA_Matches_30m, TRAFX_Biking, by = c("Location"))%>%
  dplyr::rename(TRAFx_Count = Daily_Count)

TRAFx_STRAVA_Biking_Data3 <- full_join(TRAFx_STRAVA_Biking_Data2, STRAVA_Data, by = c("date", "EdgeUID" = "edge_uid")) %>%
  dplyr::filter(!is.na(Daily_Ride_Sum))%>%
  dplyr::filter(!is.na(TRAFx_Count))

## number of TRAFx :STRAVA comparisons for biking
n_distinct(TRAFx_STRAVA_Biking_Data3$EdgeUID)

# pedestrian activities 
TRAFx_STRAVA_Pedestrian_Data2 <- full_join(TRAFx_STRAVA_Matches_30m, TRAFx_Pedestrian, by = c("Location"))%>%
  dplyr::rename(TRAFx_Count = Daily_Count)

TRAFx_STRAVA_Pedestrian_Data3 <- full_join(TRAFx_STRAVA_Pedestrian_Data2, STRAVA_Data, by = c("date", "EdgeUID" = "edge_uid")) %>%
  dplyr::filter(!is.na(Daily_Ped_Sum))%>%
  dplyr::filter(!is.na(TRAFx_Count))

## number of TRAFx :STRAVA comparisons
n_distinct(TRAFx_STRAVA_Pedestrian_Data3$EdgeUID)



# PEARSON'S CORRELATION ANALYSES ####

# Analysis 2.3.1 - Comparison of Monthly Counts Across Time  ----------------------------------------------------

#### ALL ACTIVITIES #####
##### Camera : TRAFx #####
# monthly sums, grouped by location, year and month
Cam_TRAFx_AllActivities_Monthly <- Cam_TRAFx_200m_Data3 %>%
  dplyr::group_by(Cam_Location, year, month) %>%
  dplyr::summarize(Cam_Monthly = sum (Cam_Daily), TRAFx_Monthly = sum(TRAFx_Fixed)) %>%
  dplyr::ungroup()

# calculate correlation per year/month, using scaled count values and absolute correlation
Cam_TRAFx_Monthly_Cor <- Cam_TRAFx_AllActivities_Monthly %>%
  group_by(year, month) %>%
  dplyr::summarize(cor = cor(scale(Cam_Monthly), scale(TRAFx_Monthly))) %>%
  mutate(Abs_Cor = abs(cor)) %>%
  ungroup() %>%
  dplyr::filter(!is.na(Abs_Cor))

# monthly correlation plot
## create plot of correlation by month, data grouped into year with a trend line
B1 <- ggplot(data = Cam_TRAFx_Monthly_Cor, aes(x = as.factor(month), y = Abs_Cor,color = as.factor(year), group = as.factor(year))) +
  geom_line(linewidth = 1) +
  scale_color_manual(values=c("#999999", "#5b5b5b", "#1e1e1e")) +
  geom_point(size = 3, aes(shape = as.factor(year), color = as.factor(year)))+
  geom_smooth(method = "gam", aes(group = 1),se = F, colour = 'red', linewidth = 1) +
  ylim(0.0, 1.0) +
  labs( y = "Correlation", x = "Month of the year",title = "a)") +
  theme_bw()
## for interactive plot
ggplotly(B1)


##### Camera : STRAVA #####
# monthly sums, grouped by location, year and month
Cam_STRAVA_AllActivites_Monthly <- Cam_STRAVA_Data %>%
  dplyr::group_by(Location, year, month) %>%
  dplyr::summarize(Cam_Monthly = sum (Cam_Daily), STRAVA_Monthly = sum(STRAVA_Daily)) %>%
  dplyr::ungroup()

# calculate correlation per year/month, using scaled count values and absolute correlation
Cam_STRAVA_Monthly_Cor <- Cam_STRAVA_AllActivites_Monthly %>%
  group_by(year, month) %>%
  dplyr::summarize(cor = cor(scale(Cam_Monthly), scale(STRAVA_Monthly))) %>%
  mutate(Abs_Cor = abs(cor)) %>%
  ungroup() %>%
  dplyr::filter(!is.na(Abs_Cor))

# monthly correlation plot
# create plot of correlation by month,  data grouped into year with a trend line
B3 <- ggplot(data = Cam_STRAVA_Monthly_Cor, aes(x = as.factor(month), y = Abs_Cor,color = as.factor(year), group = as.factor(year))) +
  geom_line(linewidth = 1) +
  scale_color_manual(values = c("#999999", "#5b5b5b", "#1e1e1e"))+
  geom_point(size = 3,aes(shape = as.factor(year), color = as.factor(year)))+
  geom_smooth(method = "gam", aes(group = 1),se = F, colour = 'red', linewidth = 1) +
  ylim(0.0, 1.0) +
  labs(y = "Correlation", x = "Month of the year",title = "d)") +
  theme_bw() 

# for interactive plot
ggplotly(B3)


##### TRAFx : STRAVA #####
# monthly sums, grouped by location, year and month
TRAFx_STRAVA_AllActivites_Monthly <- TRAFX_STRAVA_AllData6 %>%
  dplyr::group_by(Location, Date_Year, Date_Month) %>%
  dplyr::summarize(TRAFx_Monthly = sum (TRAFx_Fixed), STRAVA_Monthly = sum(STRAVA_Daily)) %>%
  dplyr::ungroup()

# calculate correlation per year/month, using scaled count values and absolute correlation
TRAFx_STRAVA_Monthly_Cor <- TRAFx_STRAVA_AllActivites_Monthly %>%
  group_by(Date_Year, Date_Month) %>%
  dplyr::summarize(cor = cor(scale(TRAFx_Monthly), scale(STRAVA_Monthly))) %>%
  mutate(Abs_Cor = abs(cor)) %>%
  ungroup() %>%
  dplyr::filter(!is.na(Abs_Cor)) %>%
  dplyr::filter(Date_Year < 2020)

# monthly correlation plot
# create plot of correlation by month,  data grouped into year with a trend line
B5 <- ggplot(data = TRAFx_STRAVA_Monthly_Cor, aes(x = as.factor(Date_Month), y = Abs_Cor,color = as.factor(Date_Year), group = as.factor(Date_Year))) +
  geom_line(linewidth = 1)+
  scale_color_manual(values = c("#999999", "#5b5b5b", "#1e1e1e")) +
  geom_point(size = 3,aes(shape = as.factor(Date_Year), color = as.factor(Date_Year))) +
  geom_smooth(method = "gam", aes(group = 1),se = F, colour = 'red', linewidth = 1) +
  ylim(0.0, 1.0) +
  labs(y = "Correlation", x = "Month of the year", title = "g)") +
  theme_bw() 

# for interactive plot
ggplotly(B5)


# arranging all activity plots
Plots <- ggarrange(B1 + rremove("ylab") + rremove("xlab"), B3 + rremove("ylab") + rremove("xlab"), B5 + rremove("ylab") + rremove("xlab"), ncol =1, common.legend = TRUE,legend = "none")  
Plotsb <- annotate_figure(Plots,
                          top = text_grob("All activities"))
#ggsave("Temporal_AllActivities_Monthly_Cor.jpeg", height = 5, width = 5)


### PEDESTRIAN ACTIVITIES ####
##### Camera : TRAFx #####
# monthly sums, grouped by location, year and month
Cam_TRAFx_Pedestrian_Monthly <- Cam_TRAFx_PedestrianData %>%
dplyr::group_by(Cam_Location, TRAFx_Location, year, month) %>%
  dplyr::summarize(Cam_Monthly = sum (Cam_Daily), TRAFx_Monthly = sum(TRAFx_Daily)) %>%
  dplyr::ungroup()
# calculate correlation per year/month, using scaled count values and absolute correlation
Cam_TRAFx_Ped_Monthly_Cor <- Cam_TRAFx_Pedestrian_Monthly %>%
  group_by(year, month) %>%
  dplyr::summarize(cor = cor(scale(Cam_Monthly), scale(TRAFx_Monthly))) %>%
  mutate(Abs_Cor = abs(cor)) %>%
  ungroup() %>%
  dplyr::filter(!is.na(Abs_Cor))

# monthly correlation plot
# create plot of correlation by month,  data grouped into year with a trend line
C1 <- ggplot(data = Cam_TRAFx_Ped_Monthly_Cor, aes(x = as.factor(month), y = Abs_Cor,color =as.factor(year), group = as.factor(year))) +
  geom_line(linewidth = 1) +
  scale_color_manual(values=c("#999999", "#5b5b5b", "#1e1e1e")) +
  geom_point(size = 3,aes(shape = as.factor(year), color = as.factor(year))) +
  geom_smooth(method = "gam", aes(group = 1),se = F, colour = 'red', linewidth = 1)+ 
  labs( y ="Correlation", x = "Month of the year",title = "b)") +
  ylim(0.0, 1.0) +
  theme_bw()

# for interactive plot
ggplotly(C1)


##### Camera : STRAVA #####
# monthly sums, grouped by location, year and month
Cam_STRAVA_Pedestrian_Monthly <- Cam_STRAVA_Pedestrian_Data3 %>%
  dplyr::group_by(Location, year, month) %>%
  dplyr::summarize(Cam_Monthly = sum (Cam_Daily), STRAVA_Monthly = sum(Daily_Ped_Sum)) %>%
  dplyr::ungroup() %>%
  dplyr::filter(!is.na(Cam_Monthly)) %>%
  dplyr::filter(!is.na(STRAVA_Monthly))

# calculate correlation per year/month, using scaled count values and absolute correlation
Cam_STRAVA_Ped_Monthly_Cor <- Cam_STRAVA_Pedestrian_Monthly %>%
  group_by(year, month) %>%
  dplyr::summarize(cor = cor(scale(Cam_Monthly), scale(STRAVA_Monthly))) %>%
  mutate(Abs_Cor = abs(cor)) %>%
  ungroup() %>%
  dplyr::filter(!is.na(Abs_Cor))

# monthly correlation plot
# create plot of correlation by month,  data grouped into year with a trend line
C3 <- ggplot(data = Cam_STRAVA_Ped_Monthly_Cor, aes(x = as.factor(month), y = Abs_Cor, color = as.factor(year), group = as.factor(year))) +
  geom_line(linewidth = 1)+
  scale_color_manual(values = c("#999999", "#5b5b5b", "#1e1e1e")) +
  geom_point(size = 3,aes(shape=as.factor(year), color=as.factor(year)))+
  geom_smooth(method = "loess", aes(group=1),se = F, colour = 'red', linewidth = 1) + 
  ylim(0.0, 1.0) +
  labs( y ="Correlation", x = "Month of the year",title = "e)") +
  ylim(0.0, 1.0) +
  theme_bw()

# for interactive plot
ggplotly(C3)


##### TRAFx:STRAVA #####
# monthly sums, grouped by location, year and month
TRAFx_STRAVA_Pedestrain_Monthly <- TRAFx_STRAVA_Pedestrian_Data3 %>%
  dplyr::group_by(Location, Date_Year, Date_Month) %>%
  dplyr::summarize(TRAFx_Monthly = sum(TRAFx_Count), STRAVA_Monthly = sum(Daily_Ped_Sum)) %>%
  dplyr::ungroup() %>%
  dplyr::filter(! is.na(TRAFx_Monthly)) %>%
  dplyr::filter(! is.na(STRAVA_Monthly))

# calculate correlation per year/month, using scaled count values and absolute correlation
TRAFx_STRAVA_Pedestrian_Monthly_Cor <- TRAFx_STRAVA_Pedestrain_Monthly  %>%
  group_by(Date_Year, Date_Month) %>%
  dplyr::summarize(cor=cor(scale(STRAVA_Monthly), scale(TRAFx_Monthly))) %>%
  mutate(Abs_Cor = abs(cor)) %>%
  ungroup() %>%
  dplyr::filter(!is.na(Abs_Cor)) %>%
  dplyr::filter(Date_Year < 2020)

# monthly correlation plot
# create plot of correlation by month,  data grouped into year with a trend line
C5 <- ggplot(data = TRAFx_STRAVA_Pedestrian_Monthly_Cor,  aes(x = as.factor(Date_Month), y = Abs_Cor,color = as.factor(Date_Year), group = as.factor(Date_Year))) +
  geom_line(linewidth = 1) +
  scale_color_manual(values = c("#999999", "#5b5b5b", "#1e1e1e"))+
  geom_point(size = 3,aes(shape = as.factor(Date_Year), color = as.factor(Date_Year))) +
  geom_smooth(method = "gam", aes(group = 1),se = F, colour = 'red', linewidth = 1) + 
  labs( y = "Correlation", x = "Month of the year",title = "h)") +
  ylim(0.0, 1.0) +
  theme_bw()

# cor interactive plot
ggplotly(C5)


# arranging all pedestrian plots
Plots2 <- ggarrange(C1 + rremove("ylab") + rremove("xlab"), C3 + rremove("ylab") + rremove("xlab"), C5 + rremove("ylab") + rremove("xlab"), ncol = 1, common.legend = TRUE,legend= "none")  
Plots2b <- annotate_figure(Plots2,
                           top = text_grob("Pedestrian"))
#ggsave("Temporal_Pedestrian_Monthly_Cor.jpeg", height = 5, width = 5)



### BIKING ACTIVITIES ####
#### Camera : TRAFx #####
# monthly sums, grouped by location, year and month
Cam_TRAFx_Biking_Monthly <- Cam_TRAFx_BikingData %>%
  dplyr::group_by(Cam_Location, TRAFx_Location, year, month) %>%
  dplyr::summarize(Cam_Monthly = sum(Cam_Daily), TRAFx_Monthly = sum(TRAFx_Daily), Sample_Size = n()) %>%
  dplyr::ungroup()

# calculate correlation per year/month, using scaled count values and absolute correlation
Cam_TRAFx_Biking_Monthly_Cor <- Cam_TRAFx_Biking_Monthly %>%
  group_by(year, month) %>%
  dplyr::summarize(cor = cor(scale(Cam_Monthly), scale(TRAFx_Monthly))) %>%
  mutate(Abs_Cor = abs(cor)) %>%
  ungroup() %>%
  dplyr::filter(! is.na(Abs_Cor))

# monthly correlation plot
# create plot of correlation by month,  data grouped into year with a trend line
D1 <- ggplot(data = Cam_TRAFx_Biking_Monthly_Cor,  aes(x = as.factor(month), y = Abs_Cor,color =as.factor(year), group = as.factor(year))) +
  geom_line(linewidth =1) +
  scale_color_manual(values = c("#999999", "#5b5b5b", "#1e1e1e")) +
  geom_point(size = 3, aes(shape = as.factor(year), color = as.factor(year))) +
  geom_point(size = 2, pch = 21, colour = "black",alpha = 0.7 ) +
  geom_smooth(method = "loess", aes(group = 1),se = F, colour = 'red', linewidth = 1) +
  ylim(0.0, 1.0) +
  labs( y ="Correlation", x = "Month of the year", title = "c)") +
  theme_bw()

# for interactive plot
ggplotly(D1)

#### Camera  : STRAVA #####
# monthly sums, grouped by location, year and month
Cam_STRAVA_Biking_Monthly <- Cam_STRAVA_Biking_Data3 %>%
  dplyr::group_by(Location, year, month) %>%
  dplyr::summarize(Cam_Monthly = sum (Cam_Bike_Daily), STRAVA_Monthly = sum(Daily_Ride_Sum), Sample_Size = n()) %>%
  dplyr::ungroup() %>%
  dplyr::filter(! is.na(Cam_Monthly)) %>%
  dplyr::filter(! is.na(STRAVA_Monthly))

# calculate correlation per year/month, using scaled count values and absolute correlation
Cam_STRAVA_Biking_Monthly_Cor <- Cam_STRAVA_Biking_Monthly %>%
  group_by(year, month) %>%
  dplyr::summarize(cor = cor(scale(Cam_Monthly), scale(STRAVA_Monthly)), MeanDaySample_Size = mean(Sample_Size), Monthly_sampleSize = n()) %>%
  mutate(Abs_Cor = abs(cor)) %>%
  ungroup() %>%
  dplyr::filter(! is.na(Abs_Cor))

# monthly correlation plot
# create plot of correlation by month,  data grouped into year with a trend line
D3 <- ggplot(data = Cam_STRAVA_Biking_Monthly_Cor, aes(x = as.factor(month), y = Abs_Cor,color =as.factor(year), group = as.factor(year))) +
  geom_line(linewidth = 1) +
  scale_color_manual(values = c("#999999", "#5b5b5b", "#1e1e1e")) +
  geom_point(size = 3,aes(shape=as.factor(year), color = as.factor(year))) +
  geom_smooth(method = "loess", aes(group = 1),se = F, colour = 'red', linewidth = 1) +
  ylim(0.0, 1.0) +
  labs(y ="Correlation", x = "Month of the year", title = "f)") +
  theme_bw()

# for interactive plot
ggplotly(D3)


#### TRAFx : STRAVA #####
# monthly sums, grouped by location, year and month# Monthly sums, grouped by location, year and month
TRAFx_STRAVA_Biking_Monthly <- TRAFx_STRAVA_Biking_Data3 %>%
  dplyr::group_by(Location, Date_Year, Date_Month) %>%
  dplyr::summarize(TRAFx_Monthly = sum (TRAFx_Count), STRAVA_Monthly = sum(Daily_Ride_Sum )) %>%
  dplyr::ungroup() %>%
  dplyr::filter(! is.na(TRAFx_Monthly)) %>%
  dplyr::filter(! is.na(STRAVA_Monthly))

# calculate correlation per year/month, using scaled count values and absolute correlation
TRAFx_STRAVA_Biking_Monthly_Cor <- TRAFx_STRAVA_Biking_Monthly %>%
  group_by(Date_Year, Date_Month) %>%
  dplyr::summarize(cor=cor(scale(TRAFx_Monthly), scale(STRAVA_Monthly))) %>%
  mutate(Abs_Cor = abs(cor)) %>%
  ungroup() %>%
  dplyr::filter(! is.na(Abs_Cor)) %>%
  dplyr::filter(Date_Year < 2020)

# monthly correlation plot
# create plot of correlation by month,  data grouped into year with a trend line
g <- guide_legend("Year")
D5 <- ggplot(data = TRAFx_STRAVA_Biking_Monthly_Cor, aes(x = as.factor(Date_Month), y = Abs_Cor,color = as.factor(Date_Year), group = as.factor(Date_Year))) +
  geom_line(linewidth = 1)+
  scale_color_manual(values = c("#999999", "#5b5b5a", "#1e1e1e")) +
  geom_point(size = 3,aes(shape=as.factor(Date_Year), color = as.factor(Date_Year))) +
  geom_smooth(method = "loess", aes(group = 1),se = F, colour = 'red', linewidth = 1) + 
  ylim(0.0, 1.0) +
  labs( y ="Correlation", x = "Month of the year", title = "i)") + #guides(colour = g, size = g, shape = g) +
 theme_bw() + theme(legend.position = "none") #change legend position to top, to get a legend

# for interactive plot
ggplotly(D5)

# arrangling all biking plots
Plots3 <- ggarrange(D1 + rremove("ylab") + rremove("xlab"), D3 + rremove("ylab") + rremove("xlab"), D5 + rremove("ylab") + rremove("xlab"), ncol = 1, common.legend = TRUE,legend = "none")  
Plots3b <- annotate_figure(Plots3,
                           top = "Biking")
#ggsave("Temporal_Biking_Monthly_Cor_20230814_loess.jpeg", height = 5, width = 5)

# arranging all plots
PlotsA <- ggarrange(Plotsb, Plots2b, Plots3b, nrow =1, common.legend = TRUE ) 
annotate_figure(PlotsA,
                left = text_grob("Correlation", rot = 90),
                bottom = text_grob("Month of the year"))

PlotsA_patchwork <- (plot_spacer()+ Plotsb + plot_spacer()+ Plots2b + plot_spacer() + Plots3b) +
  plot_layout(nrow = 1, widths =  c(1, 40, 1, 40, 1, 40))


zA <- wrap_elements(panel = PlotsA_patchwork) +
  labs(tag = "Correlation") +
  theme(
    plot.tag = element_text(size = rel(1.2), angle =90),
    plot.tag.position = c(0.02, 0.5)) 

Z1A <- (plot_spacer() + zA) +
  plot_layout(nrow = 1, widths =  c(0.01, 40))


Z2A <- wrap_elements(panel = Z1A) +
  labs(tag = "Month of the year") +
  theme(
    plot.tag = element_text(size = rel(1.2)),
    plot.tag.position = c(0.55, 0.03)
  )

Z1aA <- (Z2A + plot_spacer()) +
  plot_layout(nrow = 1, widths =  c(40, 0.01))

Z1aA2 <- wrap_elements(panel = Z1aA) +
  labs(tag = "Camera:Trail counter            Camera:Strava Metro         Trail counter:Strava Metro") +
  theme(
    plot.tag = element_text(size = rel(1.2), angle =270),
    plot.tag.position = c(0.97, 0.48)) 

#ggsave("Temporal_Cor_20231221.jpeg", height = 7, width = 10)
## the legend is added on top of the plot outside of R 



# Analysis 2.3.1 - Comparison of Monthly Counts Across Space -----------

### ALL ACTIVITIES ####
#### Camera : TRAFx #### 
# calculate correlation per location, using scaled count values and absolute correlation
Cam_TRAFx_AllActivities_Cor_Location <- Cam_TRAFx_AllActivities_Monthly %>%
group_by(Cam_Location) %>%
  dplyr::summarize(cor = cor(scale(Cam_Monthly), scale(TRAFx_Monthly))) %>%
  mutate(Abs_Cor = abs(cor)) %>%
  ungroup() %>%
  dplyr::filter(! is.na(Abs_Cor))

# location correlation plot
# create plot showing distribution of correlation values,  with a vertical mean line
E1a <- ggplot(data = Cam_TRAFx_AllActivities_Cor_Location, aes(x = Abs_Cor)) +
  geom_histogram(colour = "black", fill = "white", bins = 100) +
  geom_density(alpha = .2, fill = "#999999") +
  labs(x = "Correlation", y = "Count", title = "a)") +
  theme_bw() +
  theme(axis.title = element_text(size = 12), axis.text = element_text(size = 12)) +
  geom_vline(aes(xintercept = round(mean(Abs_Cor), 2)),
             color = "red", linetype = "dashed", linewidth = 1) +
  geom_text(aes(label = round(mean(Abs_Cor), digits = 2), y = 10, x = mean(Abs_Cor), vjust = -1, hjust = 1.2)) +
  scale_x_continuous(limits = c(0, NA), breaks = seq(0, 1.00, by = 0.25))

#### Camera : STRAVA #####
# calculate correlation per location, using scaled count values and absolute correlation
Cam_STRAVA_AllActivities_Cor_Location <- Cam_STRAVA_AllActivites_Monthly %>%
  group_by(Location) %>%
  dplyr::summarize(cor = cor(scale(Cam_Monthly), scale(STRAVA_Monthly))) %>%
  mutate(Abs_Cor = abs(cor)) %>%
  ungroup() %>%
  dplyr::filter(! is.na(Abs_Cor))

# location correlation plot
# create plot showing distribution of correlation values,  with a vertical mean line
E2a <- ggplot(data = Cam_STRAVA_AllActivities_Cor_Location, aes(x = Abs_Cor)) +
  geom_histogram(colour = "black", fill = "white", bins = 100) +
  geom_density(alpha = .2, fill = "#999999") +
  labs(x = "Correlation", y = "Count", title = "d)") +
  theme_bw() +
  theme(axis.title = element_text(size = 12), axis.text = element_text(size = 12)) +
  geom_vline(aes(xintercept = round(mean(Abs_Cor), 2)),
             color = "red", linetype = "dashed", linewidth = 1) +
  geom_text(aes(label = round(mean(Abs_Cor), digits = 2),y = 10,x = mean(Abs_Cor), vjust = -1, hjust = 1.2)) +
  scale_x_continuous(limits = c(0, NA), breaks = seq(0,1.00, by = 0.25))


#### TRAFx : STRAVA ####
# calculate correlation per location, using scaled count values and absolute correlation
TRAFx_STRAVA_AllActivities_Cor_Location <- TRAFx_STRAVA_AllActivites_Monthly %>%
  group_by(Location) %>%
  dplyr::summarize(cor = cor(scale(TRAFx_Monthly), scale(STRAVA_Monthly))) %>%
  mutate(Abs_Cor = abs(cor))%>%
  ungroup() %>%
  dplyr::filter(!is.na(Abs_Cor))

# location correlation plot
# create plot showing dsitribution of correlation values,  with a vertical mean line
E3a <- ggplot(data = TRAFx_STRAVA_AllActivities_Cor_Location, aes(x = Abs_Cor)) +
  geom_histogram(colour = "black", fill = "white", bins = 100) +
  geom_density(alpha = .2, fill="#999999") +
  labs(x="Correlation", y = "Count", title = "g)") +
  theme_bw() +
  theme( axis.title = element_text(size = 12),axis.text = element_text(size = 12)) +
  geom_vline(aes(xintercept = round(mean(Abs_Cor), 2)),
             color="red", linetype = "dashed", linewidth = 1) +
  geom_text(aes(label = round(mean(Abs_Cor), digits = 2), y = 10, x = mean(Abs_Cor), vjust = -1, hjust = 1.2)) +
  scale_x_continuous(limits = c(0, NA), breaks = seq(0, 1.00, by = 0.25))

# combining all plots
Plots4 <- ggarrange(E1a + rremove("ylab") + rremove("xlab"), E2a + rremove("ylab") + rremove("xlab"), E3a + rremove("ylab") + rremove("xlab"), ncol = 1)  
Plots4A <- annotate_figure(Plots4,
                           top = "All activities")
#ggsave("Spatial_AllActivities_Monthly_Cor_Histograms.jpeg", height = 6, width = 5)


## PEDESTRIAN ####
##### Camera : TRAFx ####
# calculate correlation per location, using scaled count values and absolute correlation
Cam_TRAFx_Pedestrian_Cor_Location <- Cam_TRAFx_Pedestrian_Monthly %>%
  group_by(Cam_Location) %>%
  dplyr::summarize(cor = cor(scale(Cam_Monthly), scale(TRAFx_Monthly))) %>%
  mutate(Abs_Cor = abs(cor)) %>%
  ungroup() %>%
  dplyr::filter(! is.na(Abs_Cor))

# location correlation plot
# create plot showing dsitribution of correlation values,  with a vertical mean line
F1a <- ggplot(data = Cam_TRAFx_Pedestrian_Cor_Location, aes(x = Abs_Cor)) +
  geom_histogram(colour = "black", fill = "white", bins = 100) +
  geom_density(alpha = .2, fill = "#999999") +
  labs(x = "Correlation", y = "Count", title = "b)") +
  theme_bw() +
  theme(axis.title = element_text(size = 12),axis.text = element_text(size = 12)) +
  geom_vline(aes(xintercept = round(mean(Abs_Cor), 2)),
             color="red", linetype = "dashed", linewidth = 1) +
  geom_text(aes(label = round(mean(Abs_Cor),digits = 2), y = 10,x = mean(Abs_Cor), vjust = -1, hjust = 1.2)) +
  scale_x_continuous(limits = c(0, NA),breaks = seq(0, 1.00, by = 0.25))

##### Camera : STRAVA ####
# calculate correlation per location, using scaled count values and absolute correlation
Cam_STRAVA_Pedestrian_Cor_Location <- Cam_STRAVA_Pedestrian_Monthly %>%
  group_by(Location) %>%
  dplyr::summarize(cor = cor(scale(Cam_Monthly), scale(STRAVA_Monthly))) %>%
  mutate(Abs_Cor = abs(cor)) %>%
  ungroup() %>%
  dplyr::filter(! is.na(Abs_Cor))

# location correlation plot
# create plot showing dsitribution of correlation values,  with a vertical mean line
F2a <- ggplot(data = Cam_STRAVA_Pedestrian_Cor_Location, aes(x = Abs_Cor)) +
  geom_histogram(colour = "black", fill = "white", bins = 100) +
  geom_density(alpha = .2, fill = "#999999") +
  labs(x = "Correlation", y = "Count", title = "e)") +
  theme_bw() + 
  theme(axis.title = element_text(size = 12),axis.text = element_text(size = 12)) +
  geom_vline(aes(xintercept = round(mean(Abs_Cor), 2)),
             color = "red", linetype = "dashed", linewidth = 1) +
  geom_text(aes(label = round(mean(Abs_Cor), digits = 2), y = 10, x = mean(Abs_Cor), vjust = -1, hjust = 1.2)) +
  scale_x_continuous(limits = c(0, NA),breaks = seq(0, 1.00, by = 0.25))


#### TRAFx : STRAVA ####
# calculate correlation per location, using scaled count values and absolute correlation
TRAFx_STRAVA_Pedestrian_Cor_Location <- TRAFx_STRAVA_Pedestrain_Monthly %>%
  group_by(Location) %>%
  dplyr::summarize(cor = cor(scale(TRAFx_Monthly), scale(STRAVA_Monthly))) %>%
  mutate(Abs_Cor = abs(cor)) %>%
  ungroup() %>%
  dplyr::filter(! is.na(Abs_Cor))

# location correlation plot
# create plot showing dsitribution of correlation values,  with a vertical mean line
F3a <- ggplot(data = TRAFx_STRAVA_Pedestrian_Cor_Location, aes(x = Abs_Cor)) +
  geom_histogram(colour = "black", fill = "white", bins = 100) +
  geom_density(alpha = .2, fill = "#999999") +
  labs(x = "Correlation", y = "Count", title = "h)") +
  theme_bw() +
  theme(axis.title = element_text(size = 12),axis.text = element_text(size = 12)) +
  geom_vline(aes(xintercept = round(mean(Abs_Cor), 2)),
             color = "red", linetype = "dashed", linewidth = 1) +
  geom_text(aes(label = round(mean(Abs_Cor), digits = 2), y = 6, x = mean(Abs_Cor), vjust = -1, hjust = 1.2)) +
  scale_x_continuous(limits = c(0, NA), breaks = seq(0, 1.00, by = 0.25))

# combining all plots
Plots5 <- ggarrange(F1a + rremove("ylab") + rremove("xlab"), F2a + rremove("ylab") + rremove("xlab"), F3a + rremove("ylab") + rremove("xlab"), ncol = 1, common.legend = TRUE, legend="right")  
Plots5A <- annotate_figure(Plots5,
                           top = "Pedestrian")
#ggsave("Spatial_Pedestrian_Monthly_Cor_Histograms.jpeg", height = 6, width = 5)


## BIKING Spatial ####
##### Camera : TRAFx ####
# calculate correlation per location, using scaled count values and absolute correlation
Cam_TRAFx_Biking_Cor_Location <- Cam_TRAFx_Biking_Monthly %>%
  group_by(Cam_Location) %>%
  dplyr::summarize(cor=cor(scale(Cam_Monthly), scale(TRAFx_Monthly))) %>%
  mutate(Abs_Cor = abs(cor))%>%
  ungroup() %>%
  dplyr::filter(!is.na(Abs_Cor))

# location correlation plot
# create plot showing distribution of correlation values,  with a vertical mean line
G1a <- ggplot(data = Cam_TRAFx_Biking_Cor_Location, aes(x = Abs_Cor)) +
  geom_histogram(colour = "black", fill = "white", bins = 100) +
  geom_density(alpha = .2, fill = "#999999") +
  labs(x = "Correlation", y = "Count", title = "c)") +
  theme_bw() +
  theme(axis.title = element_text(size = 12), axis.text = element_text(size = 12)) +
  geom_vline(aes(xintercept = round(mean(Abs_Cor), 2)),
             color = "red", linetype = "dashed", linewidth = 1) +
  geom_text(aes(label = round(mean(Abs_Cor), digits = 2), y = 2, x = mean(Abs_Cor), vjust = -1, hjust = 1.2)) +
  scale_x_continuous(limits = c(0, NA), breaks = seq(0, 1.00, by = 0.25))

#### Camera : STRAVA #####
# calculate correlation per location, using scaled count values and absolute correlation
Cam_STRAVA_Biking_Cor_Location <- Cam_STRAVA_Biking_Monthly %>%
  group_by(Location) %>%
  dplyr::summarize(cor = cor(scale(Cam_Monthly), scale(STRAVA_Monthly)), Sample_Size = n()) %>%
  mutate(Abs_Cor = abs(cor)) %>%
  ungroup() %>%
  dplyr::filter(! is.na(Abs_Cor))

# location correlation plot
# create plot showing distribution of correlation values,  with a vertical mean line
G2a <- ggplot(data = Cam_STRAVA_Biking_Cor_Location, aes(x = Abs_Cor)) +
  geom_histogram(colour = "black", fill = "white", bins = 100) +
  geom_density(alpha=.2, fill="#999999")+
  labs(x="Correlation", y = "Count", title = "f)")+
  theme_bw() +
  theme(axis.title = element_text(size = 12),axis.text = element_text(size = 12)) +
  geom_vline(aes(xintercept = round(mean(Abs_Cor), 2)),
             color = "red", linetype = "dashed", linewidth = 1) +
  geom_text(aes(label = round(mean(Abs_Cor), digits = 2), y = 10, x = mean(Abs_Cor), vjust = -1, hjust = 1.2)) +
  scale_x_continuous(limits = c(0, NA), breaks = seq(0, 1.00, by = 0.25))

#### TRAFx : STRAVA #####
# calculate correlation per location, using scaled count values and absolute correlation
TRAFx_STRAVA_Biking_Cor_Location <- TRAFx_STRAVA_Biking_Monthly %>%
  group_by(Location) %>%
  dplyr::summarize(cor = cor(scale(TRAFx_Monthly), scale(STRAVA_Monthly))) %>%
  mutate(Abs_Cor = abs(cor)) %>%
  ungroup() %>%
  dplyr::filter(! is.na(Abs_Cor))

# location correlation plot
# create plot showing distribution of correlation values,  with a vertical mean line
G3a <- ggplot(data = TRAFx_STRAVA_Biking_Cor_Location, aes(x = Abs_Cor)) +
  geom_histogram(colour = "black", fill = "white", bins = 100) +
  geom_density(alpha = .2, fill = "#999999") +
  labs(x = "Correlation", y = "Count", title = "i)") +
  theme_bw() +
  theme(axis.title = element_text(size = 12), axis.text = element_text(size = 12)) +
  geom_vline(aes(xintercept = round(mean(Abs_Cor), 2)),
             color = "red", linetype = "dashed", linewidth = 1) +
  geom_text(aes(label = round(mean(Abs_Cor), digits = 2), y = 3, x = mean(Abs_Cor), vjust = -1,hjust = 1.2)) +
  scale_x_continuous(limits = c(0, NA), breaks = seq(0, 1.00, by = 0.25))

# arranging all plots
Plots6 <- ggarrange(G1a + rremove("ylab") + rremove("xlab"), G2a + rremove("ylab") + rremove("xlab"), G3a + rremove("ylab") + rremove("xlab"), ncol = 1, common.legend = TRUE, legend="right")  
Plots6A <- annotate_figure(Plots6,
                           top = "Biking")
Plots7_patchwork <- (plot_spacer() + Plots4A + plot_spacer() + Plots5A + plot_spacer() + Plots6A) +
  plot_layout(nrow = 1, widths =  c(1, 40, 1, 40 ,1, 40))

# adding axes labers to combined plots
z <-wrap_elements(panel = Plots7_patchwork) +
  labs(tag = "Frequency of values") +
  theme(
    plot.tag = element_text(size = rel(1.2), angle = 90),
    plot.tag.position = c(0.02, 0.5)) 

Z1<-(plot_spacer() + z) +
  plot_layout(nrow = 1, widths =  c(0.01, 40))


Z2 <- wrap_elements(panel = Z1) +
  labs(tag = "Correlation") +
  theme(
    plot.tag = element_text(size = rel(1.2)),
    plot.tag.position = c(0.55, 0.03)
  )

Z1a <-(Z2 + plot_spacer()) +
  plot_layout(nrow = 1, widths =  c(40, 0.01))

wrap_elements(panel = Z1a) +
  labs(tag = "Camera:Trail counter            Camera:Strava Metro         Trail counter:Strava Metro") +
  theme(
    plot.tag = element_text(size = rel(1.2), angle = 270),
    plot.tag.position = c(0.97, 0.48)) 


#ggsave("PedestrianBiking_Activities_Matches20231019.jpeg", height = 7, width = 10)
# ggsave(wrap_elements(panel = Z1a) +
#          labs(tag = "Camera:Trail counter            Camera:Strava Metro         Trail counter:Strava Metro") +
#          theme(
#            plot.tag = element_text(size = rel(1.2), angle = 270),
#            plot.tag.position = c(0.97, 0.48)),
#        
#        filename = "Spatial_Corrlealtion_202314_V2.jpeg",
#        path = file.path(getwd(),"Plots"),
#        height = 8, width = 13, dpi = 300)


# Analysis 2.3.2 Comparison of Annual Counts Across Space -----------
# Comparing Strava Heatmap to annual median of trail cameras, trail counters (TRAFx) and Strava Metro counts 


# STRAVA HEATMAP GRAPHS #
## data for these graphs were generated through ArcPro, where the maximum heatmap value within 50m of a point was joined to the point through "spatial join". 

# set cut off value for number of minimum days 
## this value balances STRAVA Heatmap data representing cumulative counts over a year and including the most cameras and counters in the analysis  
cut.value <- 273 # is nine months. 

### Cameras ####
# load heatmap : camera data
Camera_HeatmapSPJoin <- read.csv("./Count_Data_MS/Camera_Heatmap_Max_SpJn_20230622.csv") %>%
  dplyr::select(Location ,UTMe, UTMn, grid_code) %>%
  dplyr::mutate(location = tolower(Location)) %>%
  dplyr::rename(easting = UTMe, northing = UTMn)

# load camera data
Camera_Data_Heatmap <- read.csv("./Count_Data_MS/CameraData_Daily_PCA_ABPARKS_CovInt.csv") # clean up columns to include all counts (not separated by activity type/species)

# remove any locations where camera was not running long enough (need at least 9 months of data - 273 days - of data x the number of years cameras were deployed)
Camera_Data_Heatmap_filter <- unique(Camera_Data_Heatmap[,c(5,7,22,24)]) %>% # get the number of active days (camera deployment) by camera location across the four years (cameras deployed: 2019, 2018, 2017, 2016) 
  dplyr::group_by(location) %>% 
  dplyr::summarise(
    # sum of total days active from the start to end date of images 
    daysactive_sum = sum(daysactive),
    # number of years between start date and last end date of image (i.e., years sampled)
    nyr_startend = length(unique(c(unique(format(as.POSIXct(lastimage, format = "%Y-%m-%d"),"%Y")), unique(format(as.POSIXct(datestart, format = "%Y-%m-%d"),"%Y")))))
  ) %>% 
  mutate(total_mindays = nyr_startend*cut.value, # for a total of 9 months of data per each year samples. 
         keepcamera = ifelse(daysactive_sum >= total_mindays, "1", "0")) %>%  # 1 if keep camera, 0 if remove 
  select(location, keepcamera) %>% 
  left_join(Camera_Data_Heatmap, by = "location", multiple = "all") %>% 
  filter(keepcamera == 1) # select only cameras where running long enough 

# select relevant columns, group by location and calculate daily total recreation counts
Camera_Data_Heatmap1 <- Camera_Data_Heatmap_filter %>%
  dplyr::select(location, dateimage, totaldaily, species, easting, northing) %>%
  group_by(location, dateimage, easting, northing) %>% 
  dplyr::summarise(TotalDaily = sum(c(totaldaily), na.rm = TRUE)) %>%
  ungroup() %>%
  separate(dateimage, c("year", "month","day"),
           remove = FALSE) %>%
  mutate(dateimage = as.POSIXct(dateimage, format = "%Y-%m-%d"))

# summarize camera data per year
# yearly totals per location
Camera_Yearly1 <- Camera_Data_Heatmap1 %>%
  dplyr::group_by(year, location, easting, northing) %>%
  dplyr::summarise(Yearly_Count = sum(TotalDaily), Sample_Size = n()) %>%
  dplyr::ungroup()

# calculate mean and median counts per location
Camera_Yearly_Summary <- Camera_Yearly1 %>%
  dplyr::group_by(location, easting, northing) %>%
  dplyr::summarise(Yearly_Mean = mean(Yearly_Count), Yearly_Median = median(Yearly_Count)) %>%
  dplyr::ungroup()

# join heatmap and camera data
Camera_Yearly_Join <- left_join(Camera_HeatmapSPJoin, Camera_Yearly_Summary, by = c("location", "easting", "northing")) %>%
dplyr::filter(! is.na(Yearly_Mean)) %>%
  dplyr::mutate(Bins = cut_interval(grid_code, length = 20))

# create graph
A <- ggplot(data = Camera_Yearly_Join, aes(x = grid_code, y = Yearly_Median))  +
  geom_point(size = 2, colour = "darkgrey",alpha = 0.3)+
  geom_point(size = 2, pch = 21, colour = "black",alpha = 0.3)+
  theme(legend.position = "none")+
  geom_smooth(colour = "red", se = FALSE, linewidth = 0.4) +
  labs(x='Strava Heatmap colour index', y = "Camera annual median", tag = "A") + 
  xlim(0, 260) +
  scale_x_continuous(breaks = seq(0,260, by = 20)) +
  scale_y_continuous(label = comma,breaks = seq(0,80000, by = 10000)) +
  theme_bw() +
  theme(legend.position = "none", panel.grid.minor = element_blank())
A

### TRAFx ####
# load Heatmap : TRAFx data
TRAFx_HeatmapSPJoin<-read.csv("./Count_Data_MS/TRAFx_Heatmap_Max_SpJn_20230622.csv") %>%
dplyr::mutate(Location = tolower(Location)) %>%
  dplyr::select(Location, grid_code, UTMn, UTMe)

# load TRAFx data counts
TRAFx_Data_HeatMap <- read.csv("./Count_Data_MS/All_TRAFx_Counts_Daily_20230525.csv") %>%
  dplyr::select(Location, Daily_Count, date) %>%
  mutate(date = as.POSIXct(date, format = "%Y-%m-%d") )%>%
  mutate(Location = tolower(Location)) %>%
  separate(date, c("year", "month","day"),
           remove = FALSE)

# calculate yearly counts
TRAFx_Yearly <- TRAFx_Data_HeatMap %>%
  dplyr::group_by(year, Location) %>%
  dplyr::arrange(date) %>% 
  mutate(DaysActive_PerYear = as.numeric(range(date, na.rm = TRUE)[2] - range(date, na.rm = TRUE)[1])) %>% # calc. number of days device operating per year 
  dplyr::ungroup() %>%
  dplyr::filter(year < 2020) %>% # remove data from after 2020
  dplyr::filter(DaysActive_PerYear < cut.value) %>% # remove locations where device is active for less than 9 months 
  dplyr::group_by(year, Location) %>% # regroup so can summary by year and TRAFx device 
  dplyr::summarise(Yearly_Count = sum(Daily_Count), Sample_Size = n()) 

# calculate mean and median per location 
TRAFx_Yearly_Avg <- TRAFx_Yearly %>%
  dplyr::group_by(Location) %>%
  dplyr::summarise(Yearly_Mean = mean(Yearly_Count), Yearly_Median = median(Yearly_Count), Sample_Size = n()) %>%
  dplyr::ungroup()

# join heatmap values with TRAFx yearly meadian
TRAFx_Heatmap_Yearly_Counts_Join <- left_join(TRAFx_Yearly_Avg, TRAFx_HeatmapSPJoin, by = c("Location")) %>%
  dplyr::filter(! is.na(grid_code)) %>%
  dplyr::mutate(Bins = cut_interval(grid_code, length = 20))

C <- ggplot(data = TRAFx_Heatmap_Yearly_Counts_Join, aes(x = grid_code, y = (Yearly_Median))) +
  geom_point(size = 2, colour = "darkgrey",alpha = 0.3)+
  geom_point(size = 2, pch = 21, colour = "black", alpha = 0.3 ) +
  theme(legend.position = "none") +
  geom_smooth(colour = "red", se = FALSE, linewidth = 0.4) +
  labs(x = 'Strava Heatmap index', y = "Counter annual median", tag = "B") + 
  scale_y_comma(label = comma) +
  theme_bw() +
  scale_x_continuous(breaks=seq(0, 260, 20)) +
  scale_y_continuous(label=comma,breaks = seq(0,508996, 100000)) + 
  expand_limits(x = 0, y = 0) +
  theme(legend.position = "none", panel.grid.minor = element_blank()) 
C


### STRAVA Metro ####
# load monthly downloads of strava data
Metro_Monthly_Counts <- read.csv("./Count_Data_MS/Monthly_Download_STRAVA_20220923.csv")

# load point locations where Strava Metro and heatmap overlapped (max heatmap value within 50m)
STRAVA_Locations <- read.csv("Count_Data_MS/Metro_WesternCA_Heatmap_MaxSPJOIN_20230621.csv") %>%
  dplyr::select(KGedgeUID  ,PreKGedgeUID , pointid ) 

# load Metro heatmap values (generated in ArcMap through spatial join to the location points)
Metro_HeatmapSPJoin  <- read.csv("./Count_Data_MS/Metro_WesternCA_Heatmap_MaxSPJOIN_20230621.csv") %>%
  dplyr::select(KGedgeUID,PreKGedgeUID, grid_code )
Metro_HeatmapSPJoin <- Metro_HeatmapSPJoin %>% 
  dplyr::filter(duplicated(Metro_HeatmapSPJoin[,c(1,2,3)]) == FALSE)

# this Strava Metro data set contains two separate sets of edge_IUDs (one from an early version of downloads, and the other from a download after an update in 2022) - some of these are for the same spatial location. One set may represent biking activities, and the other pedestrian.
# join counts to locations from one set of edge_uids
Metro_Heatmap_MonthlyDataA <-inner_join(STRAVA_Locations, Metro_Monthly_Counts , by = c("KGedgeUID" = "edge_uid")) 
# join counts to locations from another set of edge_uds
Metro_Heatmap_MonthlyDataB <-inner_join(STRAVA_Locations, Metro_Monthly_Counts , by = c("PreKGedgeUID" = "edge_uid")) 
# combine both sets
Metro_Monthly_CountsC <- bind_rows(Metro_Heatmap_MonthlyDataA, Metro_Heatmap_MonthlyDataB) 

# sum counts (bike and pedestrian) for individual pointid - these represent specific 50m radius areas that are being compared to heatmap values
Metro_Monthly_CountsD <- Metro_Monthly_CountsC%>%
  dplyr:: group_by(KGedgeUID, PreKGedgeUID, Date_Year, Date_Month, month) %>%
  dplyr:: summarise(Total_Monthly = sum(Monthly_Tot_Count)) %>%
  dplyr:: ungroup()

# calculate yearly totals
STRAVA_Yearly_Total <- Metro_Monthly_CountsD %>%
  dplyr::group_by(Date_Year, KGedgeUID, PreKGedgeUID) %>%
  dplyr::summarise(Yearly_total = sum(Total_Monthly), Sample_Size = n()) %>%
  dplyr::ungroup() %>%
  dplyr::filter(Date_Year < 2020) #remove data from after 2020
## not filter by 9 months because Strava Metro 'collects' data continuously, unlike cameras/counters with deployment periods (when active)

# calculate yearly means and medians by grouping by edgeuids
Metro_Yearly_Means_Medians <- STRAVA_Yearly_Total %>%
  dplyr::group_by(KGedgeUID, PreKGedgeUID) %>%
  dplyr::summarise(Yearly_Means = mean(Yearly_total), Yearly_Median = median(Yearly_total))%>%
  dplyr::ungroup()
# join heatmap values to median using edgeuids
Metro_Heatmap_Yearly_Mean_Join_12MonthsA <-inner_join(Metro_Yearly_Means_Medians, Metro_HeatmapSPJoin , by = c("KGedgeUID","PreKGedgeUID")) 

# remove any duplicated entries
Metro_Heatmap_Yearly_Mean_Join_12MonthsB <-Metro_Heatmap_Yearly_Mean_Join_12MonthsA %>%
  dplyr::filter(!duplicated(Metro_Heatmap_Yearly_Mean_Join_12MonthsA[,c(1,2)]))

# create graph
B <- ggplot(data = Metro_Heatmap_Yearly_Mean_Join_12MonthsB, aes(x = grid_code, y = (Yearly_Median))) +
  geom_point(size = 2, colour = "darkgrey",alpha = 0.3) +
  geom_point(size = 2, pch = 21, colour = "black",alpha = 0.3 ) +
  geom_smooth(colour = "red", se = FALSE, linewidth = 0.4) + 
  labs( x='Strava Heatmap index', y = "Strava Metro annual median", tag = "C") +
  scale_y_continuous(label=comma,breaks = seq(0,125000, 20000)) +
  theme_bw()+
  xlim(0, 260) + 
  theme(legend.position = "none", panel.grid.minor = element_blank())+
  scale_x_continuous(limits = c(0, NA),breaks = seq(0, 260, by = 20))
B

### Combine all figures together ####
PlotsA <- ggarrange(A + rremove("xlab"), C + rremove("xlab"), B + rremove("xlab"), ncol = 1) 
PlotsA <-annotate_figure(PlotsA,
                         bottom = text_grob("Strava Heatmap index", size = 11),
                         left = text_grob("Frequency of values", rot = 90))


PlotsA + 
  plot_annotation(tag_levels = list(c('A', 'B', "C"))) & 
  theme(plot.tag.position = c(0, 1),
        plot.tag = element_text(size = 10, hjust = 3, vjust = 0))

PlotsA
# ggsave(filename = "Heatmap_FrequencyOfValues.jpeg",
#        path = file.path(getwd(),"Plots"),
#        height = 7, width = 5, dpi = 400)



# Supplemental Information Figures -----------------------------------------

#### SI Figure S6 ####
# Camera : Trafx biking trend
Cam_TRAFx_Biking_Monthly$year_month <- as.yearmon(paste(Cam_TRAFx_Biking_Monthly$year, Cam_TRAFx_Biking_Monthly$month), "%Y %m")

A4a <- ggplot(data = Cam_TRAFx_Biking_Monthly, aes(x = year_month , y = Cam_Monthly)) +
  geom_point(shape = 1) +
  geom_smooth(method = "gam", se = F, colour = 'red', linewidth = 1)+
  labs( title ="Camera", x = "Date") +
  guides(x =  guide_axis(angle = 45)) +
  scale_y_continuous(label=comma, limits = c(0,1000))+
  theme_bw()+
  theme(plot.title = element_text( size = 13), legend.position = "none",axis.title = element_text(size = 12) , axis.text.y = element_text(size = 12),axis.text.x = element_text(size = 12))

A4b <- ggplot(data = Cam_TRAFx_Biking_Monthly, aes(x = year_month , y = TRAFx_Monthly)) +
  geom_point(shape = 1) +
  geom_smooth(method = "gam", se = F, colour = 'red', linewidth = 1)+
  labs( title ="Trail counter", x = "Date")+
  guides(x =  guide_axis(angle = 45)) +
  scale_y_continuous(label=comma, limits = c(0,20000))+
  theme_bw()+
  theme(plot.title = element_text( size = 13), legend.position = "none",axis.title = element_text(size = 12) , axis.text.y = element_text(size = 12),axis.text.x = element_text(size = 12))

GBike <- ggarrange(A4a + rremove("ylab") + rremove("xlab"), A4b + rremove("ylab")+ rremove("xlab"), ncol = 1,common.legend = TRUE, legend="right" )

# Camera : Strava biking trends
Cam_STRAVA_Biking_Monthly$year_month <- as.yearmon(paste(Cam_STRAVA_Biking_Monthly$year, Cam_STRAVA_Biking_Monthly$month), "%Y %m")

A6a <- ggplot(data = Cam_STRAVA_Biking_Monthly, aes(x = year_month , y = Cam_Monthly)) +
geom_point(shape = 1) +
  geom_smooth(method = "gam", se = F, colour = 'red', linewidth = 1)+
  labs(title ="Camera", x = "Date") +
  guides(x =  guide_axis(angle = 45)) +
  scale_y_continuous(label=comma, limits = c(0,2000))+
  theme_bw()+
  theme(plot.title = element_text( size = 13), legend.position = "none",axis.title = element_text(size = 12) , axis.text.y = element_text(size = 12),axis.text.x = element_text(size = 12))

A6b <- ggplot(data = Cam_STRAVA_Biking_Monthly, aes(x = year_month , y = STRAVA_Monthly)) +
  geom_point(shape = 1) +
  geom_smooth(method = "gam", se = F, colour = 'red', linewidth = 1)+
  labs( title ="Strava metro", x = "Date")+
  guides(x =  guide_axis(angle = 45)) +
  scale_y_continuous(label=comma, limits = c(0,2000))+
  theme_bw()+
  theme(plot.title = element_text( size = 13), legend.position = "none",axis.title = element_text(size = 12) , axis.text.y = element_text(size = 12),axis.text.x = element_text(size = 12))

DBike <-ggarrange(A6a + rremove("ylab") + rremove("xlab"), A6b+ rremove("ylab") + rremove("xlab"), ncol = 1,common.legend = TRUE, legend="right" )

# TRAFx : STRAVA biking trends
TRAFx_STRAVA_Biking_Monthly$year_month <- as.yearmon(paste(TRAFx_STRAVA_Biking_Monthly$Date_Year, TRAFx_STRAVA_Biking_Monthly$Date_Month), "%Y %m")

A9a <- ggplot(data = TRAFx_STRAVA_Biking_Monthly, aes(x = year_month , y = TRAFx_Monthly)) +
  geom_point(shape = 1) +
  geom_smooth(method = "gam", se = F, colour = 'red', linewidth = 1)+
  labs( title ="Trail counter", x = "Date")+
  guides(x =  guide_axis(angle = 45)) +
  scale_y_continuous(label=comma, limits = c(0,15000))+
  theme_bw()+
  theme(plot.title = element_text( size = 13), legend.position = "none",axis.title = element_text(size = 12) , axis.text.y = element_text(size = 12),axis.text.x = element_text(size = 12))

A9b <- ggplot(data = TRAFx_STRAVA_Biking_Monthly, aes(x = year_month , y = STRAVA_Monthly)) +
  geom_point(shape = 1) +
  geom_smooth(method = "gam", se = F, colour = 'red', linewidth = 1)+
  labs( title ="Strava metro", x = "Date")+
  guides(x =  guide_axis(angle = 45)) +
  scale_y_continuous(label=comma, limits = c(0,1000))+
  theme_bw()+
  theme(plot.title = element_text( size = 13), legend.position = "none",axis.title = element_text(size = 12) , axis.text.y = element_text(size = 12),axis.text.x = element_text(size = 12))

CBike <-ggarrange(A9a+ rremove("ylab") + rremove("xlab"), A9b+ rremove("ylab") + rremove("xlab"), ncol = 1,common.legend = TRUE, legend="right" )

# Camera : TRAFx pedestrian trends
Cam_TRAFx_Pedestrian_Monthly$year_month <- as.yearmon(paste(Cam_TRAFx_Pedestrian_Monthly$year, Cam_TRAFx_Pedestrian_Monthly$month), "%Y %m")

A3a <-ggplot(data = Cam_TRAFx_Pedestrian_Monthly, aes(x = year_month , y = Cam_Monthly)) +
  geom_point(shape = 1) +
  geom_smooth(method = "gam", se = F, colour = 'red', linewidth = 1)+
  labs( title ="Camera", x = "Date") +
  guides(x =  guide_axis(angle = 45)) +
  scale_y_continuous(label=comma, limits = c(0,15000))+
  theme_bw()+
  theme(plot.title = element_text( size = 13), legend.position = "none",axis.title = element_text(size = 12) , axis.text.y = element_text(size = 12),axis.text.x = element_text(size = 12))

A3b <-ggplot(data = Cam_TRAFx_Pedestrian_Monthly, aes(x = year_month , y = TRAFx_Monthly)) +
  geom_point(shape = 1) +
  geom_smooth(method = "gam", se = F, colour = 'red', linewidth = 1)+
  labs(title ="Trail counter", x = "Date") +
  guides(x =  guide_axis(angle = 45)) +
  scale_y_continuous(labels=comma, limits = c(0,15000))+
  theme_bw()+
  theme(plot.title = element_text( size = 13), legend.position = "none",axis.title = element_text(size = 12) , axis.text.y = element_text(size = 12),axis.text.x = element_text(size = 12))

G <- ggarrange(A3a+ rremove("ylab")+ rremove("xlab"), A3b+ rremove("ylab") + rremove("xlab"), ncol = 1,common.legend = TRUE, legend="right" )

# Camera : STRAVA pedestrian trends
Cam_STRAVA_Pedestrian_Monthly$year_month <- as.yearmon(paste(Cam_STRAVA_Pedestrian_Monthly$year, Cam_STRAVA_Pedestrian_Monthly$month), "%Y %m")

A7a <-ggplot(data = Cam_STRAVA_Pedestrian_Monthly, aes(x = year_month , y = Cam_Monthly)) +
  geom_point(shape = 1) +
  geom_smooth(method = "gam", se = F, colour = 'red', linewidth = 1)+
  labs(title ="Camera", x = "Date") +
  guides(x =  guide_axis(angle = 45)) +
  scale_y_continuous(label=comma, limits = c(0,30000))+
  theme_bw()+
  theme(plot.title = element_text( size = 13), legend.position = "none",axis.title = element_text(size = 12) , axis.text.y = element_text(size = 12),axis.text.x = element_text(size = 12))

A7b <-ggplot(data = Cam_STRAVA_Pedestrian_Monthly, aes(x = year_month , y = STRAVA_Monthly)) +
  geom_point(shape = 1) +
  geom_smooth(method = "gam", se = F, colour = 'red', linewidth = 1)+
  labs( title ="Strava metro", x = "Date")+
  guides(x =  guide_axis(angle = 45)) +
  scale_y_continuous(label=comma, limits = c(0,1000))+
  theme_bw()+
  theme(plot.title = element_text( size = 13), legend.position = "none",axis.title = element_text(size = 12) , axis.text.y = element_text(size = 12),axis.text.x = element_text(size = 12))

D <- ggarrange(A7a+ rremove("ylab") + rremove("xlab"), A7b+ rremove("ylab") + rremove("xlab"), ncol = 1,common.legend = TRUE, legend="right" )

# TRAFx : STRAVA pedestrian trends
TRAFx_STRAVA_Pedestrain_Monthly$year_month <- as.yearmon(paste(TRAFx_STRAVA_Pedestrain_Monthly$Date_Year , TRAFx_STRAVA_Pedestrain_Monthly$Date_Month ), "%Y %m")

A10a <-ggplot(data = TRAFx_STRAVA_Pedestrain_Monthly, aes(x = year_month , y = TRAFx_Monthly)) +
geom_point(shape = 1) +
  ylim(0,20000)+
  geom_smooth(method = "gam", se = F, colour = 'red', linewidth = 1)+
  labs( title ="Trail counter", x = "Date")+
  guides(x =  guide_axis(angle = 45)) +
  scale_y_continuous(label=comma)+
  theme_bw()+
  theme(plot.title = element_text( size = 13), legend.position = "none",axis.title = element_text(size = 12) , axis.text.y = element_text(size = 12),axis.text.x = element_text(size = 12))

A10b <-ggplot(data = TRAFx_STRAVA_Pedestrain_Monthly, aes(x = year_month , y = STRAVA_Monthly)) +
  geom_point(shape = 1) +
  geom_smooth(method = "gam", se = F, colour = 'red', linewidth = 1)+
  labs(title ="Strava metro", x = "Date")+
  guides(x =  guide_axis(angle = 45)) +
  scale_y_continuous(label=comma, limits=c(1,5000))+
  theme_bw()+
  theme(plot.title = element_text( size = 13), legend.position = "none",axis.title = element_text(size = 12) , axis.text.y = element_text(size = 12),axis.text.x = element_text(size = 12))

C <- ggarrange(A10a+ rremove("ylab") + rremove("xlab"), A10b+ rremove("ylab") + rremove("xlab"), ncol = 1,common.legend = TRUE, legend="right" )

# combining all into one figure
p_lab <-ggplot() + 
  annotate(geom = "text", x = 1, y = 1, label = "Total counts per location", angle = 90) +
  coord_cartesian(clip = "off")+
  theme_void()

Pedestrian_patchwork <- (plot_spacer()+G+plot_spacer()+D+plot_spacer()+C)+
  plot_layout(nrow = 1, widths =  c(1,40,1,40,1,40))

Pedestrian_patchwork2 <- Pedestrian_patchwork + 
  plot_annotation(title = "Pedestrian",
                  tag_levels = list(c('A', 'B', 'C')))  & 
  theme(plot.tag.position = c(0, 1),
        plot.tag = element_text(size = 12, hjust = 1, vjust = 0))


Pedestrian_patchwork3 <- wrap_elements(panel = Pedestrian_patchwork2) +
  labs(tag = "Total counts per loction") +
  theme(
    plot.tag = element_text(size = rel(1), angle =90),
    plot.tag.position = c(0.02, 0.5)) 


Bike_patchwork <- (plot_spacer()+GBike+plot_spacer()+DBike+plot_spacer()+CBike)+
  plot_layout(nrow = 1, widths =  c(1,40,1,40,1,40))

Bike_patchwork2 <- Bike_patchwork +  
  plot_annotation(title = "Biking" ,
                  tag_levels = list(c('D', 'E', 'F')))  & 
  theme(plot.tag.position = c(0, 1),
        plot.tag = element_text(size = 12, hjust = 1, vjust = 0))

Bike_patchwork3 <- wrap_elements(panel = Bike_patchwork2) +
  labs(tag = "Total counts per loction") +
  theme(
    plot.tag = element_text(size = rel(1), angle =90),
    plot.tag.position = c(0.02, 0.5)) 


patchwork <-(plot_spacer()/Pedestrian_patchwork3 /plot_spacer() / Bike_patchwork3)+
  plot_layout(ncol = 1, heights =  c(1,15,1,15))


all_patchwork <- wrap_elements(patchwork) + 
  labs(tag = "Date") +
  theme(
    plot.tag = element_text(size = rel(1)),
    plot.tag.position = c(0.5, 0.02)
  ) 

# png(filename = file.path(getwd(),"Plots", "FigureS6.png"),
#     width = 4000, height = 4800, res = 400)
# plot(all_patchwork)
# dev.off()


# End Script #