nztamap.R

# site locations
# https://opendata-nzta.opendata.arcgis.com/datasets/b90f8908910f44a493c6501c3565ed2d_0/explore?location=-40.549980%2C173.382464%2C6
# download of full data as csv as of 2026-04-02, approx 385KB


# traffic counters
# https://opendata-nzta.opendata.arcgis.com/datasets/898e77d403c443c4961ada073d144735_0/explore
# download of full data as csv as of 2026-04-01, approx 800MB

# road outline map 
# https://data.linz.govt.nz/layer/50329-nz-road-centrelines-topo-150k/
# full download, kml, 257.5 MB

# coastlines
# https://data.linz.govt.nz/layer/51153-nz-coastlines-and-islands-polygons-topo-150k/
# full download, kml, 20.3 MB

# rivers
# https://data.linz.govt.nz/layer/50328-nz-river-polygons-topo-150k/
# full download, kml, 25 MB

# lakes
# https://data.linz.govt.nz/layer/50293-nz-lake-polygons-topo-150k/
# full download, kml, 4 MB

library(readr)
library(dplyr)
library(lubridate)
library(ggplot2)
library(sf)
library(ggrepel)
source("~/theme.R")
library(patchwork)

#day before the weeks of interest start
time26p1t1 = ymd("2026-02-21")
time26p2t1 = ymd("2026-03-28")

#calculated times
time26p1t2 = time26p1t1 + days(8)
time25p1t1 = time26p1t1 - days(365)
time25p1t2 = time25p1t1 + days(8)
time26p2t2 = time26p2t1 + days(8)
time25p2t1 = time26p2t1 - days(365)
time25p2t2 = time25p2t1 + days(8)


# TMS locations already in NZTM coordinates
TMS_locations <- read_csv("../data/State_highway_traffic_monitoring_sites.csv")

# syncronising 2025 readings with 2026 by days of the week
TMS_readings <- read_csv("../data/TMS_Telemetry_Sites_-1825614806134071233.csv",
                         col_types=cols(
                           OBJECTID = col_double(),
                           `Start Date` = col_datetime(format = "%m/%d/%Y %I:%M:%S %p"),
                           `Site Alias` = col_double(),
                           `Region Name` = col_character(),
                           `Site Reference` = col_character(),
                           `Class Weight` = col_character(),
                           `Site Description` = col_character(),
                           `Lane Number` = col_double(),
                           `Flow Direction` = col_double(),
                           `Traffic Count` = col_double()
                         )) |> 
  mutate(matchday = if_else(year(`Start Date`) == 2026, `Start Date`, `Start Date` - days(1)))


coast <- read_sf("../data/lds-nz-coastlines-and-islands-polygons-topo-150k-KML/nz-coastlines-and-islands-polygons-topo-150k.kml")
coastline_nztm <- st_transform(coast, 2193)

roads <- read_sf("../data/lds-nz-road-centrelines-topo-150k-KML/nz-road-centrelines-topo-150k.kml")
roads_nztm <- st_transform(roads, 2193)

rivers <- read_sf("../data/lds-nz-river-polygons-topo-150k-KML/nz-river-polygons-topo-150k.kml")
rivers_nztm <- st_transform(rivers, 2193)

lakes <- read_sf("../data/lds-nz-lake-polygons-topo-1250k-KML/nz-lake-polygons-topo-1250k.kml")
lakes_nztm <- st_transform(lakes, 2193)

useful_sites <- TMS_readings |> 
  filter((matchday > time26p1t1 & matchday < time26p1t2) |
           (matchday > time26p2t1 & matchday < time26p2t2) |
           (matchday > time25p1t1 & matchday < time25p1t2) |
           (matchday > time25p2t1 & matchday < time25p2t2) ,
         `Class Weight` == "Light") |> 
  count(`Site Reference`, matchday) |> 
  count(`Site Reference`) |> 
  filter(n==28) |> 
  select(siteref = `Site Reference`)


quantdata <- TMS_readings |> 
  filter(`Site Reference` %in% useful_sites$siteref,
         (matchday > time26p1t1 & matchday < time26p1t2) |
           (matchday > time26p2t1 & matchday < time26p2t2) |
           (matchday > time25p1t1 & matchday < time25p1t2) |
           (matchday > time25p2t1 & matchday < time25p2t2) ,
         `Class Weight` == "Light") |> 
  summarise(.by=c(`Site Reference`,matchday),
            totcounts = sum(`Traffic Count`)) |> 
  arrange(`Site Reference`,matchday) |> 
  summarise(.by=`Site Reference`,
            t25a = mean(totcounts[1:7]),
            t25b=mean(totcounts[8:14]),
            t26a = mean(totcounts[15:21]),
            t26b=mean(totcounts[22:28]),
            calc = t26b/t26a - t25b/t25a,
            change= paste0(round(100 * calc,1),"%")) |> 
  rename(siteref=`Site Reference`) |> 
  inner_join(TMS_locations, by = join_by(siteref))
  

#bounding lats & longsx
longmin = 170.40074  # min longitude
latmin = -45.90207 # min latitude
longmax = 170.60153  # max longitude
latmax = -45.82076   # max latitude

# easier to check the end map bounds online in lat/long so set it up and convert
bbox_ll <- st_bbox(
  c(
    xmin = longmin,  # min longitude
    ymin = latmin,  # min latitude
    xmax = longmax,  # max longitude
    ymax = latmax   # max latitude
  ),
  crs = 4326
)
bbox_nztm <- st_bbox(st_transform(st_as_sfc(bbox_ll), 2193))

# clipping the map data makes plotting faster, particularly when experimenting 
# with the plot look, and stops the repel calcs going mad
coastc <- st_crop(coastline_nztm, bbox_nztm)
roadc <- st_crop(roads_nztm, bbox_nztm)
riverc <-  st_crop(rivers_nztm, bbox_nztm)
lakec <-  st_crop(lakes_nztm, bbox_nztm)
quantc <- quantdata |> 
  filter(X >= bbox_nztm["xmin"], X <= bbox_nztm["xmax"],
         Y >= bbox_nztm["ymin"], Y <= bbox_nztm["ymax"])

#map bounds as water underlay
bbox_geom <- st_as_sfc(bbox_nztm)

watercolour = "slategray1"
m1 <- ggplot() +
  geom_sf(data = bbox_geom, fill = watercolour, colour = NA, linewidth = 1)+
  geom_sf(data=coastc, aes(geometry=geometry),  fill="white", colour="white", lwd=.05) +
  geom_sf(data=riverc, aes(geometry=geometry),  fill=watercolour, colour=watercolour, lwd=.1) +
  geom_sf(data=lakec, aes(geometry=geometry),  fill=watercolour, colour=watercolour, lwd=.1) +
  geom_sf(data=roadc, aes(geometry=geometry),  fill=NA, colour="black", lwd=.04) +
  geom_point(data=quantc, aes(x=X,y=Y), size=0.5, colour="red") +
  coord_sf(xlim = c(bbox_nztm["xmin"], bbox_nztm["xmax"]),
           ylim = c(bbox_nztm["ymin"], bbox_nztm["ymax"])) + 
  theme_david_map() +
  theme(plot.subtitle = element_text(size=6))+
  geom_label_repel(data=quantc, aes(x=X,y=Y, label=change),
                   max.overlaps = Inf, size=2, fill="#FFFFFF99",
                   point.padding = 0.1, box.padding = .5,
                   min.segment.length = 0) +
  labs(subtitle=paste0("100*((week to ", time26p2t2 - days(1), "/week to 2026-2-28)-\n(same change for pervious year))"))

#trend
trendata <- TMS_readings |> 
  filter(`Site Reference` %in% quantc$siteref,
         (matchday >= ymd_h("2026-02-01 00") & matchday < time26p2t2) |
           (matchday >= ymd_h("2025-02-01 00") & matchday < time25p2t2) ,
         `Class Weight` == "Light") |> 
  summarise(.by = c(`Site Reference`,matchday), dcount=sum(`Traffic Count`)) |> 
  mutate(Year = year(matchday)) |> 
  arrange(Year,`Site Reference`, matchday) |> 
  group_by(Year, `Site Reference`) %>%
  mutate(Value_loess = loess(dcount ~ as.numeric(matchday))$fitted) %>%
  ungroup() |> 
  summarise(.by=c(Year, matchday),
            Trend = mean(Value_loess)) |> 
  mutate(Year = factor(year(matchday)),
         Dateis = ISOdate(2026,month(matchday), day(matchday)))

g1 <- ggplot(trendata, aes(x=Dateis, y=Trend, colour=Year)) +
  geom_line() + theme_david() +
  coord_cartesian(ylim=c(0,NA)) +
  scale_colour_manual(values=six_cols[2:1]) +
  geom_vline(xintercept = ymd("2026-2-28"), colour=six_cols[3]) +
  annotate("text", x=ymd("2026-2-27"), y=5000, label="Start of War",
           colour=six_cols[3], angle=90,vjust=0) +
  theme(legend.position = "inside",
        legend.position.inside = c(.7,.4),
        plot.subtitle = element_text(size=6)) +
  labs(subtitle="Mean Traffic Counts trend, mean(loess(site daily totals))",
       y=NULL, x=NULL)

chart <- m1 + g1 +
  plot_annotation(title="Change in light vehicle traffic in Dunedin area",
                  caption=make_footer("Sources: NZTA open data, LINZ base maps"))

ggsave(filename="~/Desktop/ggsky.jpg",plot = chart,
  width = 2016,
  height = 1134,
  units = "px",
)