This is the code used to generate the graph.

# scrape data
library(rvest)
library(tidyverse)
library(lubridate)
read_html("https://www.manchester.ac.uk/coronavirus/cases/") %>%
  html_nodes("table") %>% 
  html_table() %>%
  bind_rows() -> tbls

# import the early non-tabulated data
tibble(Date = seq(ymd('2020-09-21'), ymd('2020-09-26'), by='days'),
       Students = c(1, 1, 2, 5, 6, 6), 
       Staff = rep(NA, 6)) %>%
  mutate(Date = ymd(Date)) -> tblsOLD

# separate and label 2020 data
tbls %>%
  slice_head(n = 96) %>%
  mutate(Date = str_c(Date, " 2020")) %>%
  mutate(Date = dmy(Date)) -> tbls2020

# separate and label 2021 data
tbls %>%
  .[97:dim(.)[1],] %>%
  mutate(Date = str_c(Date, " 2021")) %>%
  mutate(Date = dmy(Date)) -> tbls2021

# join them all
bind_rows(list(tblsOLD, tbls2020, tbls2021)) %>%
  arrange(Date) -> df

# tidy-up
rm(list= ls()[! (ls() %in% c('df'))])

# pivot longer 
df %>%
  pivot_longer(c(Students, Staff), names_to = "pop") %>%
  mutate(pop = factor(pop, levels = c("Students", "Staff"))) %>%
  arrange(pop, Date) -> df

# compute rolling mean data
df_roll <- df %>% 
  group_by(Date) %>%
  summarise(sum = sum(value)) %>%
  mutate(value = round(as.vector(stats::filter(sum, rep(1/7, 7), sides = 2))))

# import significant dates for annotation
sigDates <- tibble(Date = parse_date(c("21/09/2020", "28/09/2020", 
                                        "05/10/2020", "07/10/2020", 
                                        "03/12/2020", "09/12/2020",
                                        "18/12/2020", "04/01/2021"),
                                      format = "%d/%m/%Y"),
                   labels = c("1st year induction starts", "wider induction starts", 
                              "teaching starts", "DfE Tier 3 starts", 
                              "travel window starts", "travel window ends",
                              "closure starts", "closure ends"))

# create plot
ggplot() +
  geom_area(data = df, mapping = aes(x = Date, y = value, fill = pop), alpha = 0.6) +
  geom_line(data = df_roll, mapping = aes(x = Date, y = value, col = "black")) +
  labs(fill = "", 
       subtitle = paste0("Data retrieved automatically from \"manchester.ac.uk/coronavirus/cases\" on ", Sys.Date()),
       x = "Date",
       caption = paste0("Mean no. of daily reported cases over the last 7 days = ", 
                        round(mean(tail(df_roll, n = 7)$sum), 1), 
                        " (~", 
                        round(sum(tail(df_roll, n = 7)$sum*(100000 / (12800 + 40250))), 1), 
                        " new cases/100,000/week)", "\n",
                        "Case rate estimates based on data from \"manchester.ac.uk/discover/facts-figures\"")) +
  ggtitle("University of Manchester COVID-19 Daily Reported Cases") +
  #theme_classic() +
  scale_x_date(breaks = df$Date[wday(df$Date,label = TRUE) == "Mon" & day(df$Date) <= 7], 
               date_labels = "%d %b",
               minor_breaks = df$Date[which(wday(df$Date)==2)],
               sec.axis = sec_axis(trans = ~ . + 0,
                                   breaks = sigDates$Date,
                                   labels = sigDates$labels)
  ) +
  scale_y_continuous( "Cases [n/day]", 
                      sec.axis = sec_axis(~ . * (100000/(12800+40250)), name = "Estimated daily case rate [n/100,000/day]") 
  ) +
  theme(axis.text.x.top = element_text(angle = 20, 
                                       hjust = 0)
  ) +
  geom_vline(xintercept = sigDates$Date, 
             linetype = "dotted") +
  scale_color_identity(name = "", 
                       guide = "legend", 
                       labels = "7-day rolling mean")