Herding Data: The Udderly Fascinating Trends in New Zealand Livestock

1 Introduction

In New Zealand, discussing livestock numbers is practically a national pastime, whether you’re interested in the agri-economy, environmental health, or the ever-popular sheep-to-people ratio.¹

In this post, I will present visualisations aimed at enhancing our understanding of the number and distribution of various livestock over time and across different regions of New Zealand.

2 Data sources

StatsNZ publishes livestock numbers as part of their Agricultural Production Survey (APS) derived data series. It’s important to note that within this single dataset, there are varying levels of spatial resolution and time series:

Total livestock numbers across NZ (1971-2019)
Regional breakdown (1990-2019)
APS hex grid breakdown (1994-2017)

Show the code

library(tidyverse)
library(sf)
library(mapgl)
library(koordinatr)
library(zoo)
library(gganimate)
library(scales)
library(isaacr)
library(pals)
library(magick)
library(billboarder)

options(scipen = 999)

# Load the data
aps17 <- st_read("data/livestock-numbers-grid-aps-2017/livestock-numbers-grid-aps-2017.shp")
aps94 <- st_read("data/livestock-numbers-grid-aps-1994/livestock-numbers-grid-aps-1994.shp")

aps17$sheepdens[aps17$grid_id == 724] <- 0
aps17$sheepdens[aps17$grid_id == 742] <- 0

livestock_numbers_raw <- get_table_as_tibble(Sys.getenv("mfe_api_key"), "mfe", "105406") |> select(-gml_id)

# Read in the shapefile of regions
regions_sf <- st_read("data/statsnz-regional-council-2023-clipped-generalised-SHP/regional-council-2023-clipped-generalised.shp") |> 
  mutate(REGC2023_2 = str_replace(REGC2023_2, " Region$", "")) |> 
  st_simplify(dTolerance = 1000, preserveTopology = FALSE)

3 Trends in Total Cattle Numbers

To start, let’s examine how total cattle numbers in New Zealand have changed over time.

Total cattle numbers peaked in 2014 with 10.4 million cattle.
Around the year 2000, the number of dairy cattle surpassed the number of beef cattle for the first time.
Dairy cattle also peaked in 2014 at 6.7 million, while beef cattle peaked much earlier in 1975 at 6.3 million and have been steadily declining since.

Show the code

cattle_numbers <- livestock_numbers_raw |>
  filter(geography_name == "New Zealand") |>
  filter(animal %in% c("Beef cattle", "Dairy cattle", "Total cattle")) |>
  group_by(animal) |>
  mutate(count = na.spline(count, na.rm = FALSE)) # Interpolate missing values via cubic spline interpolation

# Create the plot
p <- ggplot(cattle_numbers, aes(x = year, y = count, colour = animal, group = animal)) +
  geom_line() + 
  geom_point(size = 2) +
  geom_text(aes(label = animal), vjust = -0.5, hjust = 0, size = 4) +
  labs(title = "Total cattle numbers in New Zealand",
       x = "",
       y = "Number of cattle") + 
  theme_minimal() +
  expand_limits(x = c(1971, 2019 + 6), y = 0) + # Expand limits on x-axis
  theme(legend.position = "none") +
  scale_colour_brewer(palette = "Set1") +
  scale_y_continuous(breaks = seq(0, max(cattle_numbers$count, na.rm = TRUE), by = 1e6),
                     labels = scales::label_number(scale = 1e-6, suffix = "M")) + # Format y-axis labels 
  transition_reveal(year)

# Render the animation
animate(p,
        nframes = 100,
        fps = 10,
        end_pause = 30)

Tip

The graph below is interactive. You can hover over the time series to see the exact values, you can hover over the legend to highlight each animal type, and you can click on the legend to show/hide certain animals.

Show the code

# Convert count to millions for better readability in the plot
cattle_numbers2 <- cattle_numbers |>
  mutate(count_million = count / 1e6)

# Create the plots
billboarder() |>
  bb_linechart(
    data = cattle_numbers2,
    mapping = bbaes(x = year, y = count_million, group = animal)
  ) |>
  bb_legend(position = "right") |>
  bb_x_axis(tick = list(values = seq(1971, 2019, by = 5))) |>
  bb_y_axis(tick = list(format = suffix("M"))) |>
  bb_x_grid(show = TRUE) |> 
  bb_y_grid(show = TRUE) |> 
  bb_title(text = "Total Cattle Numbers in New Zealand") |>
  bb_colors_manual("Beef cattle" = "#DA0719", "Dairy cattle" = "#2D6BAA", "Total cattle" = "#40A33B") |>
  bb_legend(show = TRUE)

4 Decline in Sheep Numbers

Sheep numbers in New Zealand have experienced precipitous decline.

Since their peak in 1982 at just over 70 million, sheep numbers have been declining almost every year, reaching a record low of 26.8 million in 2019.

Show the code

# Load and prepare your data
sheep_numbers <- livestock_numbers_raw |>
  filter(geography_name == "New Zealand") |>
  filter(animal == "Sheep") |>
  group_by(animal) |>
  mutate(count = na.spline(count, na.rm = FALSE)) # Interpolate missing values via cubic spline interpolation

# Create the plot
p2 <- ggplot(sheep_numbers, aes(x = year, y = count, colour = animal, group = animal)) +
  geom_line() +
  geom_point(size = 2) +
  geom_text(aes(label = animal), vjust = -0.5, hjust = 0, size = 4) +
  labs(title = "Total sheep numbers for New Zealand",
       x = "",
       y = "Number of sheep") +
  theme_minimal() +
  scale_colour_brewer(palette = "Set1") +
  expand_limits(x = c(1971, 2019 + 6), y = 0) + # Expand limits on x-axis
  scale_y_continuous(breaks = seq(0, max(sheep_numbers$count, na.rm = TRUE), by = 1e7),
                     labels = scales::label_number(scale = 1e-6, suffix = "M")) + # Format y-axis labels
  theme(legend.position = "none") +
  transition_reveal(year)

# Render the animation
animate(p2, nframes = 100, fps = 10, end_pause = 20)

5 Comparing Livestock Populations

A graph of total livestock numbers highlights the sheer dominance of sheep compared to other livestock types. However, this dominance can obscure the trends in other animals, making it difficult to discern their patterns and changes over time.

Show the code

# Load and prepare your data
livestock_numbers <- livestock_numbers_raw  |> 
  filter(geography_name == "New Zealand") |>
  group_by(animal) |>
  mutate(count = na.spline(count, na.rm = FALSE)) |> # Interpolate missing values via cubic spline interpolation
  filter(count > 0) # Filter out animals with zero counts

# Create the plot
p3 <- ggplot(livestock_numbers, aes(x = year, y = count, colour = animal, group = animal)) +
  geom_line() +
  geom_point(size = 2) +
  geom_text(aes(label = animal), vjust = -0.5, hjust = 0, size = 4) +
  labs(title = "Total livestock numbers for New Zealand",
       x = "",
       y = "Number of animals") +
  theme_minimal() +
  scale_colour_brewer(palette = "Set1") +
  expand_limits(x = c(1971, 2019 + 6), y = 0) + # Expand limits on x-axis
  scale_y_continuous(breaks = seq(0, max(livestock_numbers$count, na.rm = TRUE), by = 1e7),
                     labels = scales::label_number(scale = 1e-6, suffix = "M")) + # Format y-axis labels
  theme(legend.position = "none") +
  transition_reveal(year)

# Render the animation
animate(p3, nframes = 100, fps = 10, end_pause = 20)

6 Stock Unit Equivalents Explained

The previous graph doesn’t effectively communicate the varying metabolic demand of different animals, nor the changes in metabolic demand over time. These demands have been increasing due to rising animal size and productivity.

Table 1: Stock unit equivalent values assumed for sheep, beef, dairy and deer livestock classes between 1980 and 2017. From Snelder et al (2021)

Stock type	1980	1985	2000	2017
Sheep	0.95	0.95	1.15	1.35
Beef	5	5.3	6	6.9
Dairy	5	5.5	6.8	8
Deer	1.6	1.6	2	2.3

Converting to stock unit equivalents, as done by Snelder et al. (2021),² allows for more meaningful comparisons. I’ve used their methods of conversion for consistency. In short, the count of animals is multiplied by the stock unit equivalent (Table 1) for that animal type and nearest year.

Total stock units peaked in 2012 at 121.9 million stock unit equivalents (SU).
Total stock units in New Zealand have remained relatively high over time, with a modest increase from 106 million SUs in 1980 to 115 million SUs in 2019.

Show the code

# Apply the function to the dataset
total_stock_units <- livestock_numbers_raw |> 
  filter(geography_name == "New Zealand") |>
  filter(animal != "Total cattle") |> 
  group_by(animal) |> 
  mutate(count = na.spline(count, na.rm = FALSE)) |> # Interpolate missing values via cubic spline interpolation
  ungroup() |> 
  rowwise() |>
  mutate(stock_unit_equivalent = count * get_stock_unit(animal, year)) |>
  ungroup() |> 
  group_by(year) |> 
  summarise(total_stock_unit_equivalent = sum(stock_unit_equivalent, na.rm = TRUE))

# Create the plot
p4 <- ggplot(total_stock_units, aes(x = year, y = total_stock_unit_equivalent)) +
  geom_line() +
  geom_point(size = 2) +
  geom_text(label = "Total stock units", vjust = -0.5, hjust = 0, size = 4) +
  labs(title = "Total stock unit equivalents for New Zealand",
       x = "",
       y = "Number of stock units") +
  theme_minimal() +
  expand_limits(x = c(1971, 2019 + 6), y = 0) + # Expand limits on x-axis
  scale_y_continuous(breaks = seq(0, max(total_stock_units$total_stock_unit_equivalent, na.rm = TRUE), by = 5e7),
                     labels = scales::label_number(scale = 1e-6, suffix = "M")) + # Format y-axis labels
  theme(legend.position = "none") +
  transition_reveal(year)

# Render the animation
animate(p4, nframes = 100, fps = 10, end_pause = 20)

7 Regional Livestock Trends

Although Waikato still leads significantly in terms of dairy cattle numbers, Canterbury has seen the largest increase in dairy cattle since 1990.
At their peak in 2014, Waikato had 1.9 million dairy cattle, which slightly decreased to 1.8 million by 2019.
Despite this decrease, Waikato’s dairy cattle numbers are still almost 1.5 times greater than Canterbury’s peak of 1.3 million in 2014, and 2.6 times greater than Southland’s peak of 730,000 in 2015.

Show the code

dairy_numbers_region <- livestock_numbers_raw |>
  filter(geography_type == "Region") |>
  filter(geography_name != "Chatham Islands") |> 
  filter(animal == "Dairy cattle") |>
  group_by(geography_name) |>
  mutate(count = na.approx(count, rule = 2, na.rm = FALSE)) |>  # Interpolate missing values via linear interpolation
  mutate(baseline_1990 = first(count[year == 1990])) |> 
  mutate(change_from_1990 = count - baseline_1990) |> 
  filter(count > 0) |>  # Filter out animals with zero counts |> 
  mutate(year = as.integer(year))

# Create the plot
p5 <- ggplot(dairy_numbers_region, aes(x = year, y = count, colour = geography_name, group = geography_name)) +
  geom_line() +
  geom_point(size = 2) +
  geom_text(aes(label = geography_name), vjust = -0.5, hjust = 0, size = 4) +
  labs(title = "Dairy cattle numbers by region",
       x = "",
       y = "Number of dairy cattle") +
  theme_minimal() +
  expand_limits(x = c(1990, 2019 + 6), y = 0) + # Expand limits on x-axis
  theme(legend.position = "none") +
  scale_y_continuous(labels = label_comma()) + # Format y-axis labels with commas
  scale_color_manual(values = c("Canterbury" = "red")) +
  transition_reveal(year)

# Render the animation
animate(p5,
        nframes = 100,
        fps = 10,
        end_pause = 30)

Tip

The graph below is interactive. It looks pretty busy, but you can hover over the legend to isolate individual regions. And you can click on the legend items to hide or show multiple regions.

Show the code

# Create the plots
billboarder(width = 700, height = 500) |>
  bb_linechart(
    data = dairy_numbers_region,
    mapping = bbaes(x = year, y = count, group = geography_name)
  ) |>
  bb_legend(position = "right") |>
  bb_x_axis(tick = list(values = seq(1971, 2019, by = 5))) |>
  bb_y_axis(tick = list(format = htmlwidgets::JS("function(d) { return d.toLocaleString(); }"))) |>
  bb_x_grid(show = TRUE) |> 
  bb_y_grid(show = TRUE) |> 
  bb_title(text = "Dairy cattle numbers by region") |>
  bb_legend(show = TRUE) |>
  bb_tooltip(grouped = FALSE)

8 Animated Regional Trends Map

Want the interactive version of this? Click through to this Shiny app.

As mentioned above, Canterbury stands out for its significant increase in dairy cows since 1990. In the 24 years leading up to 2014, the region added over 1.2 million dairy cattle.
At the opposite end of the spectrum is Auckland, which, by 2009, had lost 135,000 dairy cattle since 1990, reaching its lowest point.

Show the code

# Merge spatial data with dairy numbers
map_data <- dairy_numbers_region |>
  left_join(regions_sf, by = c("geography_name" = "REGC2023_2")) |> 
  filter(!is.na(LAND_AREA_)) |>  # Filter out regions with missing values |>
  #filter(geography_name %in% c("Canterbury", "Southland", "Northland", "Waikato", "Bay of Plenty")) |> 
  ungroup()

# Calculate the range of change_from_1990
dairy_change_range <- range(dairy_numbers_region$change_from_1990, na.rm = TRUE)
max_abs_change <- max(abs(dairy_change_range))

# Define the limits to be symmetric around zero
symmetric_limits <- c(-max_abs_change, max_abs_change)

# Create the line plot
line_plot <- ggplot(dairy_numbers_region, aes(x = year, y = change_from_1990, colour = change_from_1990, group = geography_name)) +
  geom_line() +
  geom_point(size = 2) +
  geom_text(aes(label = geography_name), vjust = -0.5, hjust = 0, size = 4) +
  labs(title = "Change in dairy cattle numbers from 1990, by region",
       x = "",
       y = "Number of dairy cattle") +
  theme_minimal() +
  expand_limits(x = c(1990, 2019 + 6)) + # Expand limits on y-axis to match symmetric limits
  theme(legend.position = "none") +
  scale_color_distiller(
    palette = "RdBu",
    limits = symmetric_limits,
    labels = label_comma()) +
  scale_y_continuous(labels = label_comma()) + # Format y-axis labels with commas
  transition_reveal(year)

# Create the map plot animation
map_plot <- ggplot(map_data) +
  geom_sf(aes(fill = change_from_1990, geometry = geometry), color = "black") +
  scale_fill_distiller(
    palette = "RdBu",
    limits = symmetric_limits,
    labels = label_comma(), 
    name = "Dairy cattle change \nfrom 1990") +
  theme_void() +
  theme(legend.position = "right") +
  labs(title = 'Year: {frame_time}') +
  transition_time(year) +
  ease_aes('linear')

# Render the individual animations
line_anim <- animate(line_plot, nframes = 100, fps = 10, end_pause = 30)
map_anim <- animate(map_plot, nframes = 100, fps = 10, end_pause = 30)

# Combine the animations side by side using magick
line_gif <- image_read(line_anim)
map_gif <- image_read(map_anim)
combined_gif <- image_append(c(line_gif[1], map_gif[1]))

# Ensure both GIFs have the same number of frames
num_frames_line <- nrow(image_info(line_gif))
num_frames_map <- nrow(image_info(map_gif))
num_frames <- min(num_frames_line, num_frames_map)

# Combine the animations side by side
combined_gif <- image_append(c(line_gif[1], map_gif[1]), stack = FALSE)

for (i in 2:num_frames) {
  combined <- image_append(c(line_gif[i], map_gif[i]), stack = FALSE)
  combined_gif <- c(combined_gif, combined)
}

# Add a pause by duplicating the last frame multiple times
pause_duration <- 30 # Number of frames to pause
last_frame <- combined_gif[num_frames]

for (i in 1:pause_duration) {
  combined_gif <- c(combined_gif, last_frame)
}

# Display the combined animation
combined_gif

9 2D Livestock Density Comparison

Now you get a bonus for reading this far! The talented Kyle Walker has been developing a new R package called mapgl, which brings the power of Mapbox GL JS and MapLibre GL JS to R users. This package introduces exciting new features that surpass traditional mapping packages like Leaflet.

The first map is a simple 2D representation of the livestock hex grids, allowing you to slide the slider back and forth to compare sheep density between the years 1994 and 2017.

Show the code

my_column <- "sheepdens"
my_range <- round(range(aps17$sheepdens), 0)

mapboxgl() |>
  fit_bounds(aps94, animate = FALSE) |>
  add_fill_layer(id = "APS 1994",
                 source = aps94,
                 fill_color = interpolate(
                   column = my_column,
                   values = c(0, 420), 
                   stops = c("#F4F4F3", "#A40D23"),
                   na_color = "lightgrey",
                 ),
                 fill_opacity = 0.8) |> 
  add_legend("Sheep density (count/ha) <br>1994 / 2017",
    values = my_range,
    colors = c("#F4F4F3", "#A40D23")
  )-> map1

mapboxgl() |> 
  fit_bounds(aps17, animate = FALSE) |>
  add_fill_layer(id = "APS 2017",
                 source = aps17,
                 fill_color = interpolate(
                   column = my_column,
                   values = c(0, 420), 
                   stops = c("#F4F4F3", "#A40D23"),
                   na_color = "lightgrey",
                 ),
                 fill_opacity = 0.8) -> map2

# very hacky way to include the map in the html. See issue at: https://github.com/walkerke/mapgl/issues/3
htmlwidgets::saveWidget(compare(map1, map2), "www/map1.html")

10 3D Livestock Density Visualisation

The next plot is something truly special! The hex grids are rendered in 3D, with both their height and colour proportional to the sheep density. You can slide the slider back and forth as before, but you can also hold down the control button while clicking to tilt and pan the 3D view.

In 3D, the dramatic reduction in sheep numbers is striking, and you can clearly see where the losses have predominantly occurred in the South Island.

Show the code

my_column <- "sheepdens"
my_range <- round(range(aps17$sheepdens), 0)

maplibre(
  style = maptiler_style('basic'),
  center = c(170.5, -43.5),
  zoom = 5,
  pitch = 60,
  bearing = 0
) |>
  add_fill_extrusion_layer(
    id = "3d",
    source = aps94,
    fill_extrusion_color = interpolate(
      column = my_column,
      values = my_range,
      stops = c("#F4F4F3", "#A40D23")
    ),
    fill_extrusion_height = interpolate(
      column = my_column,
      values = my_range, # Data values for interpolation
      stops = c(0, 160000) # Corresponding heights for those values
    )
  ) |>
  add_legend("Sheep density (count/ha) <br>1994 / 2017",
    values = my_range,
    colors = c("#F4F4F3", "#A40D23")
  ) -> map3

maplibre(
  style = maptiler_style('basic'),
  center = c(170.5, -43.5),
  zoom = 5,
  pitch = 60,
  bearing = 0
) |>
  add_fill_extrusion_layer(
    id = "3d",
    source = aps17,
    fill_extrusion_color = interpolate(
      column = my_column,
      values = my_range,
      stops = c('#F4F4F3', '#A40D23')
    ),
    fill_extrusion_height = interpolate(
      column = my_column,
      values = my_range,  # Data values for interpolation
      stops = c(0, 160000)  # Corresponding heights for those values
    )
  ) -> map4

#compare(map3, map4)

# very hacky way to include the map in the html. See issue at: https://github.com/walkerke/mapgl/issues/3
htmlwidgets::saveWidget(compare(map3, map4), "www/map2.html")

11 Conclusions

That wraps it up. This post explored livestock numbers in New Zealand, highlighting trends and regional variations for cattle and sheep populations. By using stock unit equivalents, we gained deeper insights into these changes.

We also showcased innovative visualisation techniques, from animations to interactive 3D comparisons, to illustrate these trends effectively.

I hope this has provided valuable insights into the livestock dataset and inspired new ways to visualise agricultural data. Thank you for reading, and stay tuned for more insights and visualisations.

Footnotes

FYI it’s currently 4.6 sheep per person! See RNZ for more details.↩︎
Snelder TH, Fraser C, Larned ST, Monaghan R, De Malmanche S, Whitehead AL. Attribution of river water-quality trends to agricultural land use and climate variability in New Zealand. Marine and Freshwater Research. 2021;73(1):1-19. doi:10.1071/mf21086↩︎