Skip to content

Instantly share code, notes, and snippets.

@thoughtfulbloke
Created March 18, 2026 03:24
Show Gist options
  • Select an option

  • Save thoughtfulbloke/e059cecd2a276e351a64644c2cddce0a to your computer and use it in GitHub Desktop.

Select an option

Save thoughtfulbloke/e059cecd2a276e351a64644c2cddce0a to your computer and use it in GitHub Desktop.
# get a clipped map of territorial authorities, like
# https://koordinates.com/from/datafinder.stats.govt.nz/layer/123496-territorial-authority-2026-clipped/
# download (needs registration) in a format useful you. I used
# Coordinate Reference System WGS 84 (EPSG:4326)) and as a .kml file (21MB)
#
# download territorial population like from
# https://infoshare.stats.govt.nz Population : Population Estimates (DPE)
# Estimated Resident Population for Territorial Authority Areas, at 30 June(1996+) (Annual-Jun)
# as a csv for the year 2025
#
# Download another dataset with similar geographies
# like https://www.stats.govt.nz/large-datasets/csv-files-for-download/
# Industries - Agricultural production statistics: key tables from 2002–2017
# agricultural-production-statistics-june-2017-final-additional-tables.xlsx
#
# put them all in the working directory for easy access
# set R to use that working directory
# load some use helper libraries
# these packages need to have been installed in a one off way first.
library(readr)
library(tidyr)
library(readxl)
library(sf)
library(dplyr)
library(ggplot2)
# read in the humans
humans <- readr::read_csv("DPE389801_20260318_023730_82.csv", skip=1, n_max = 1)
# because the csv is by default wide (1 row per year), make it long (1 row per Territory)
longform <- tidyr::pivot_longer(humans, -`...1`, names_to = "Territory", values_to = "Population")
# read in Agricultural Stats
agr <- readxl::read_excel(path = "Copy of agricultural-production-statistics-june-2017-final-additional-tables.xlsx",
sheet = "Table 2", skip=4)
# pick the relevant agriculture column
# and make the names the same as the people data
# the agricultural data has S for suppressed but there are a few so setting to 3
agrPop <- agr |> dplyr::select(Description, `TOTAL beef cattle`, `TOTAL dairy cattle`) |>
dplyr::mutate(Territory = gsub("[0123456789][0123456789][0123456789]\\. ", "", Description),
Beef = as.numeric(`TOTAL beef cattle`),
Beef = ifelse(is.na(Beef),3,Beef),
Dairy = as.numeric(`TOTAL dairy cattle`),
Dairy = ifelse(is.na(Dairy),3,Dairy)) |>
dplyr::left_join(longform, by = join_by(Territory)) |>
dplyr::filter(!is.na(Population)) #removing rows with no match
# read the map
terland <- sf::read_sf("territorial-authority-2026-clipped.kml")
# do some maths on the numeric columns to create the fill values
# change Territory names to match map name:
# Otorohanga to Ōtorohanga
# Opotiki to Ōpōtiki
# join together with the map,
# excluding chathams to solve dateline size of map problems
combo <- agrPop |>
dplyr::mutate(legs = 2*Population + 4*Beef + 4*Dairy,
Territory= gsub("Otorohanga", "Ōtorohanga", Territory),
Territory= gsub("Opotiki", "Ōpōtiki", Territory)) |>
rename(Name = Territory) |> # change data Territory to match map name
dplyr::left_join(terland, by = join_by(Name)) |>
filter(Name != "Chatham Islands Territory")
# plot the map
ggplot(combo, aes(geometry=geometry, fill=legs)) +
geom_sf(colour=NA) +
scale_fill_viridis_c()+
theme_void() +
labs(title="Human and cow legs, total, by Territorial Authroity")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment