library(plyr)  
library(tidyverse)
library(lubridate)
library(zoo)
library(scales)
library(rgdal)

# start with some labeling.
# Long names can be cumbersome during data manipulation, but their essential for when we label graphs

reg_labs <- c("Fraser", "Interior", "Northern", "Vancouver Coastal", "Vancouver Island")

# assign colours to the regions
reg_colours <- rev(RColorBrewer::brewer.pal(5,"Dark2"))
names(reg_colours) <- c("FH","INH","NH","VCH","VIH")

# make a tibble to store region names and their abbreviations

map_regions <- 
  tibble(reg2 = c("FH", "INH", "NH", "VCH", "VIH"),
         HA_Name = c("Fraser", "Interior", "Northern", "Vancouver Coastal", "Vancouver Island"))

# read in the shapefile
health_regs <- readOGR(dsn = ".", layer = "HA_2018", verbose = FALSE) 
# dsn = "." tells the function to look here in our local directory (assuming you saved your shapefile in your local directory)
# layer = "HA_2018" is the name of the shapefile, without the file extension
# verbose = false because I don't really care about seeing a progress bar as I've wandered away for another cup of coffee

# the result is a SpatialPolygonDataFrame, which we can't plot with ggplot, we need to convert this to a data frame

# we start by adding an extra id column to the SPDF
health_regs@data$id <- rownames(health_regs@data)

# then we convert the SPDF to a dataframe using fortify
# we tell fortify to use the id column that we created as our regions (this is essentially a grouping variable)
health.points <- fortify(health_regs, region = "id")

# we now have a dataframe that has all of the gps coordinates that were in the SPDF, but none of the metadata (region name, etc)

head(health.points)

# to get that info, we use the id column we created to merge the new dataframe with the metadata part of the SPDF
health_df <- join(health.points, health_regs@data, by = "id")

# add region abbreviations to this new dataframe
health_df <- 
health_df %>%
  left_join(map_regions)


# when plotting maps, we need to tell ggplot how to scale the map correctly otherwise it will just be plotted as a square
# we could tell ggplot how we want the map to be projected using a model like mercator etc. but this is good enough for what I want to do
map_ratio <- (max(health_df$lat) - min(health_df$lat))/(max(health_df$long)-min(health_df$long))

ggplot()+
  geom_polygon(data = health_df, aes(x = long, y = lat, group = group, fill = reg2), colour = "black") +
  coord_fixed(ratio = map_ratio) +
  scale_fill_manual(name = "Health Region",
                    values = reg_colours,
                    labels = reg_labs)+
  theme(axis.title = element_blank(),
        axis.text = element_blank(),
        legend.position = c(0.9, 0.9),
        axis.ticks = element_blank(),
        panel.background = element_blank())

daily_reports <- read_csv("BC_reports.csv", trim_ws = TRUE, na=c("","NA"))

daily_reports

# select just the regional data and convert to long form

reg_reports <- 
daily_reports %>%
  select(date, contains("H", ignore.case = FALSE)) %>%
  gather(key = region, value = cume_cases, -date) 

reg_reports

# store the most recent date of the data we have

most_recent <- 
reg_reports %>%
  select(date) %>%
  pull() %>%
  max(.) %>%
  as_date()

# store the length of time between the earliest date for which we don't have data, and the most recent date

fill_dates_length <- 1 + interval(as_date("2020-03-13"),most_recent) %/% ddays(1)

# create a dataframe with the full complement of dates and regions and join with our reported data
# this creates the rows where we are missing data
reg_reports_full <- 
tibble(date = rep(seq(as_date("2020-03-13"), most_recent, by = "day"),5),
       region = rep(c("FH","INH","NH","VCH","VIH"), each = fill_dates_length)) %>%
  left_join(reg_reports)

reg_reports_full

# we can now identify the gaps in the data 
gaps <- 
reg_reports_full %>%
  mutate(rownum = 1:NROW(reg_reports_full)) %>%
  filter(is.na(cume_cases)) %>%
  pull(rownum)

# for each row in the full data set that is missing case numbers, store that row number and the one before and after  
gap_brack <- vector()
for(i in 1:length(gaps)){
  gap_brack <- c(gap_brack, gaps[i]-1, gaps[i], gaps[i] + 1)
}

# keep the unique values (remove extras due to multiple rows missing data)
gap_brack <- unique(gap_brack)

# now we can select the rows that we need to interpolate for, and perform the interpolation
reg_reports_gap <- 
reg_reports_full %>%
  mutate(rownum = 1:NROW(reg_reports_full)) %>%
  filter(rownum %in% gap_brack) %>%
  mutate(cume_cases = na.approx(cume_cases)) %>%
  select(-rownum)

reg_reports_gap

# in order for this to plot well, we need to add the dates for which we do have data to the interpolated data tibble
# but we don't want the real data in this tibble, because if it is included, it will appear in the plot as interpolated data
# we create the full complement of dates and regions again, and merge our interpolated data onto it
reg_reports_gap <- 
tibble(date = rep(seq(as_date("2020-03-13"), most_recent, by = "day"),5),
       region = rep(c("FH","INH","NH","VCH","VIH"), each = fill_dates_length)) %>%
  left_join(reg_reports_gap)

reg_reports_gap

# bind the reported and interpolated data together, labelling them as such with a new column (type)
reg_reports_all <- 
bind_rows(reg_reports %>%
            add_column(type = "repo"),
          reg_reports_gap %>%
            add_column(type = "inter"))

# plot
ggplot(data = reg_reports_all, aes(x = date, y = cume_cases))+
  geom_line(aes(colour = region, linetype = type, size = type))+ # set linetype based on data type
  scale_linetype_manual(name = "Data type", values = c("dotted", "solid"), labels = c("Interpolated", "Reported"))+
  scale_size_manual(name = "Data type", values = c(0.5, 1), labels = c("Interpolated", "Reported"))+
  scale_colour_manual(name = "Health Region", values = reg_colours, labels = reg_labs)+ # set colours using our predifined vectors
  xlab("Date")+
  ylab("Cumulative # of Cases")+
  ylim(0,710)+
  theme(legend.position = c(0.15, 0.6))

reg_details <- read_csv("reg_details.csv", trim_ws = TRUE, na = c("","NA"))

reg_details

# the categories column has long strings, which is inconvenient to work with when manipulating data
# I've created abbreviatons for these categories in a separate csv

cat_labs <- read_csv("cat_labs.csv", trim_ws = TRUE, na = c("","NA"))

cat_labs

# add the labels 
reg_details <- 
reg_details %>%
  left_join(cat_labs) 

# calculate the number of active cases by subtracting cumulative recoveries and deaths from cumulative cases

reg_actives <-
reg_details %>%
  filter(cat == "cume_cases" | cat == "cume_recov" | cat == "cume_died") %>% # filter just the data we want
  select(-total, -range, -n, -comments, -category) %>% # select just the columns we need
  gather(key = reg, value = n, -date, -cat) %>% # convert to long form
  spread(key = cat, value = n) %>% # convert back to wide form, this time with case types as columns
  mutate(active = cume_cases - cume_recov - cume_died) %>% # calculate actives
  left_join(map_regions, by = c("reg" = "HA_Name")) # add region abbreviations

reg_actives

# store the most recent date
most_recent2 <- 
reg_details %>%
  select(date) %>%
  pull() %>%
  max(.) %>%
  as_date()

# store the length of the interval between the earliest and latest dates in the set
fill_dates_length2 <- 1 + interval(as_date("2020-03-23"),most_recent2) %/% ddays(1)

# create the full complement of dates and regions to create the rows that we're missing
# join with the existing date
reg_reports_full2 <- 
tibble(date = rep(seq(as_date("2020-03-23"), most_recent2, by = "day"),5),
       reg2 = rep(c("FH","INH","NH","VCH","VIH"), each = fill_dates_length2)) %>%
  left_join(reg_actives)

reg_reports_full2

# most regions have the same gaps, but Vancouver Coastal, for whatever reason, is missing more data 
# to compensate, we can do the interpolation with the other regions first

# remove Vancouver coastal
# identify gaps
gaps2 <- 
reg_reports_full2 %>% 
  filter(reg2 != "VCH") %>%
  mutate(rownum = 1:NROW(.)) %>%
  filter(is.na(cume_cases)) %>%
  pull(rownum)
  
# for each gap in the data (row with missing data), store that row and the one before and after
gap_brack2 <- vector()
for(i in 1:length(gaps2)){
  gap_brack2 <- c(gap_brack2, gaps2[i]-1, gaps2[i], gaps2[i] + 1)
}

# store just unique values
gap_brack2 <- unique(gap_brack2)

# perform the interpolation
reg_reports_gap2 <- 
reg_reports_full2 %>%
  filter(reg2 != "VCH") %>%
  mutate(rownum = 1:NROW(.)) %>%
  filter(rownum %in% gap_brack2) %>%
  mutate(cume_cases = na.approx(cume_cases),
         cume_recov = na.approx(cume_recov),
         active = na.approx(active),
         cume_died = na.approx(cume_died)) %>%
  select(-rownum) %>%
  left_join(map_regions, by = c("reg2" = "reg2")) %>%
  select(-reg) %>%
  rename(., reg = HA_Name)

# Now do vancouver coastal

gaps3 <- 
reg_reports_full2 %>% 
  filter(reg2 == "VCH") %>%
  mutate(rownum = 1:NROW(.)) %>%
  filter(is.na(active)) %>%
  pull(rownum)
  
gap_brack3 <- vector()
for(i in 1:length(gaps3)){
  gap_brack3 <- c(gap_brack3, gaps3[i]-1, gaps3[i], gaps3[i] + 1)
}

gap_brack3 <- unique(gap_brack3)

reg_reports_gap3 <- 
reg_reports_full2 %>%
  filter(reg2 == "VCH") %>%
  mutate(rownum = 1:NROW(.)) %>%
  filter(rownum %in% gap_brack3) %>%
  mutate(cume_cases = na.approx(cume_cases),
         cume_recov = na.approx(cume_recov),
         active = na.approx(active),
         cume_died = na.approx(cume_died)) %>%
  select(-rownum) %>%
  left_join(map_regions, by = c("reg2" = "reg2")) %>%
  select(-reg) %>%
  rename(reg = HA_Name)

#combine 

reg_reports_gap2 <- 
  bind_rows(reg_reports_gap2,
            reg_reports_gap3)

# add in the empty dates

reg_reports_gap2 <- 
tibble(date = rep(seq(as_date("2020-03-23"), most_recent2, by = "day"),5),
       reg2 = rep(c("FH","INH","NH","VCH","VIH"), each = fill_dates_length2)) %>%
  left_join(reg_reports_gap2)

# put the reported data and interpolated data together in a single data frame

reg_reports_all2 <- 
bind_rows(reg_reports_full2 %>%
            left_join(map_regions, by = c("reg2" = "reg2")) %>%
            select(-reg) %>%
            rename(reg = HA_Name) %>%
            add_column(type = "repo"),
          reg_reports_gap2 %>%
            left_join(map_regions, by = c("reg2" = "reg2")) %>%
            select(-reg) %>%
            rename(reg = HA_Name) %>%
            add_column(type = "inter"))

reg_reports_all2

# so here, we could use the care labels that we imported at the beginning (cat_labs)
# but I want even shorter labels than what is in the original data
# so I add my own labels, manually, and order them as a factor so that the figure legend is ordered the way I want it to be 
# ggplot2 will otherwise order legend entries alphabetically
rra2_long <- 
reg_reports_all2 %>%
  gather(key = cases, value = ntot, -date, -reg2, -reg, -type) %>%
  mutate(case_type = factor(ifelse(cases == "cume_died", "Cumulative Deaths", 
                                   ifelse(cases == "cume_recov", "Cumulative Recoveries", 
                                          ifelse(cases == "active", "Active Cases", "Cumulative Cases"))), 
                            levels = c("Cumulative Cases", "Cumulative Recoveries", "Cumulative Deaths", "Active Cases"), 
                            ordered = TRUE))

rra2_long

# there are two points in the data set that I want to highlight with asterisks, cumulative recoveries in Vancouver coastal on 2020-03-26 and 2020-04-17 

marker <- 
tibble(date = c(ymd("2020-03-26"),ymd("2020-04-17")),
       ntot = c(175,460),
       case_type = c("Cumulative Recoveries","Cumulative Recoveries")) %>%
         mutate(case_type = factor(case_type, levels = c("Cumulative Cases", "Cumulative Recoveries", "Cumulative Deaths", "Active Cases"), ordered = TRUE))


# Lastly, there are a couple spots in the data where we have a single day of data and then gaps on either side
# when we plot the reported data as a line graph, these single days will not appear
# because they are isolated and there are no adjacent points to draw lines with
# we can plot these isolated days as points as a work-around
# I identify these single points by looking at the gap data
# I calculate the difference between dates between each consequtive pair of rows
# if the difference is equal to two days (ie there is a day with data between them)
# I select the later gap date and convert it back to the date with data
# then use left join to select the correct data
max_date <-
  rra2_long %>%
  select(date) %>%
  pull() %>%
  unique() %>%
  max(.)


extra_dates <-
  bind_rows(rra2_long %>%
              filter(type == "repo" & is.na(ntot)) %>% # filter data to return just reported NAs (gaps)
              mutate(date2 = c(NA,diff(date) == 2)) %>% 
              filter(date2 == TRUE) %>%
              select(date, reg2, cases, type) %>%
              mutate(date = date - 1) %>%
              left_join(rra2_long),
            rra2_long %>%
              filter(date == max_date)) # I also add on the data from the latest date in the dataset, in case those data are also isolated

# I store the earliest date to use when plotting
min_date <-
rra2_long %>%
  select(date) %>%
  pull() %>%
  unique() %>%
  min(.)

ggplot(data = rra2_long, aes(x = date, y = ntot))+
  geom_line(aes(colour = reg2, linetype = type, size = type))+ # plot the reported and interpolated data as lines
  geom_point(data = extra_dates, aes(x = date, y = ntot, colour = reg2), size = 1, shape = 15, show.legend = FALSE)+ # add the single points
  geom_point(data = marker, aes(x = date, y = ntot), size = 2, shape = 8)+ # add the asterisks
  scale_linetype_manual(name = "Data type", values = c("dotted", "solid"), labels = c("Interpolated", "Reported"))+
  scale_size_manual(name = "Data type", values = c(0.5, 1), labels = c("Interpolated", "Reported"))+
  scale_colour_manual(name = "Health Region", values = reg_colours, labels = reg_labs)+
  xlab("Date")+
  ylab("# of Cases")+
  xlim(c(min_date, max_date))+
  facet_wrap(.~case_type, ncol = 2) +
  ylim(0,710)+
  theme(legend.position = c(0.15, 0.3))

reg_pops <- read_csv("reg_pops.csv", trim_ws = TRUE, na = c("","NA"))

reg_pops

# we may did into population structure later, for now we just want totals

reg_tots <- 
reg_pops %>%
  slice(1) %>%
  select(-cat) %>%
  gather(key = region, value = pop_tot)

ggplot(data = reg_tots, aes(x = region, y = pop_tot))+
  geom_bar(stat = "identity", aes(fill = region), width = 0.4) +
  scale_x_discrete(name = "Health Region", labels = reg_labs) +
  scale_y_continuous(name = "Population (# of People)", labels = comma) +
  scale_fill_manual(values = reg_colours, guide = FALSE)

per_cap_actives <- # save plot 
reg_reports_all2 %>%
  left_join(reg_tots, by = c("reg2" = "region")) %>% # left join the population date onto our dataframe
  mutate(per_cap = 100000*active/pop_tot) %>% # calculate cases per 100,000
  ggplot(aes(x = date, y = per_cap))+ # plot
  geom_line(aes(colour = reg2, linetype = type, size = type))+
  scale_linetype_manual(name = "Data type", values = c("dotted", "solid"), labels = c("Interpolated", "Reported"))+
  scale_size_manual(name = "Data type", values = c(0.5, 1), labels = c("Interpolated", "Reported"))+
  scale_colour_manual(name = "Health Region", values = reg_colours, labels = reg_labs)+
  xlab("Date")+
  ylim(0,20)+
  ylab("Cumulative # of Cases per 100,000 People")

# repeat the process for our single data points
extra_dates_actives <-
extra_dates %>%
  left_join(reg_tots, by = c("reg2" = "region")) %>%
  mutate(per_cap = 100000*ntot/pop_tot) %>%
  filter(cases == "active")

per_cap_actives + # add the single data points to the plot
  geom_point(data = extra_dates_actives, aes(x = date, y = per_cap, colour = reg2), 
             size = 1, shape = 15, show.legend = FALSE)