Ontario Wildfire Analysis

library(sf)
library(spdep)
library(raster)
library(sp)
library(rasterVis)
library(spatstat)
library(dbscan)
library(ggplot2)
library(ggspatial)
library(dplyr)
library(readr)
library(viridis)
library(knitr)
library(leaflet)
sf_use_s2(FALSE)

Fire_Disturbance_Point <- read.csv("Fire_Disturbance_Point.csv")
fire_data <- read_csv("Fire_Disturbance_Point.csv", show_col_types = FALSE)

# Fire cause codes as labelled in the MNRF Fire Disturbance Point data
cause_labels <- c(
  IDF = "Industrial (forest operations)", IDO = "Industrial (other)", INC = "Incendiary",
  LTG = "Lightning", MIS = "Miscellaneous", REC = "Recreation", RES = "Resident",
  RWY = "Railway", UNK = "Unknown"
)
cause_colours <- c(
  IDF = "#98971a", IDO = "#928374", INC = "#d62221", LTG = "#d79921", MIS = "#b16286",
  REC = "#458592", RES = "#8f3f71", RWY = "#de8d32", UNK = "#665c54"
)

# Basemap layers (Natural Earth, public domain), prepared by data/prepare_basemap.R
ontario     <- st_read("data/ontario.gpkg", quiet = TRUE)
great_lakes <- st_read("data/great_lakes.gpkg", quiet = TRUE)
neighbours  <- st_read("data/neighbours.gpkg", quiet = TRUE)

fires_sf <- st_as_sf(Fire_Disturbance_Point, coords = c("X", "Y"), crs = 4269, remove = FALSE)

# Ontario MNR Lambert, so the province keeps its true shape
ON_CRS <- 3161
ontario_bbox <- st_bbox(st_transform(ontario, ON_CRS))

base_map <- function() {
  list(
    geom_sf(data = neighbours, fill = "#8c857714", colour = "#8c857755", linewidth = 0.2, inherit.aes = FALSE),
    geom_sf(data = ontario, fill = "#8c857722", colour = "#8c8577", linewidth = 0.3, inherit.aes = FALSE),
    geom_sf(data = great_lakes, fill = "#45859233", colour = "#458592aa", linewidth = 0.2, inherit.aes = FALSE)
  )
}
map_frame <- function() {
  list(
    coord_sf(crs = ON_CRS, xlim = ontario_bbox[c("xmin", "xmax")], ylim = ontario_bbox[c("ymin", "ymax")],
             expand = TRUE),
    annotation_north_arrow(location = "tr", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"),
                           style = north_arrow_minimal(line_col = ms_text, fill = ms_text, text_col = ms_text)),
    annotation_scale(location = "bl", width_hint = 0.25, bar_cols = c(ms_text, "transparent"),
                     line_col = ms_text, text_col = ms_text),
    theme_mikael(),
    theme(axis.title = element_blank())
  )
}

Overview

This analysis looks at the Ontario Ministry of Natural Resources and Forestry (MNRF) Fire Disturbance Point dataset: 63,767 recorded fires from 1976 to 2022, each with a location, year, cause and final size in hectares.

The sections below cover how large fires get, how size and frequency change over the years, where fires happen, what causes them, and where large lightning-caused fires cluster.

Fire Size

Fire sizes are extremely skewed: most fires are tiny, while a handful burn tens of thousands of hectares.

size_summary <- summary(Fire_Disturbance_Point$FIRE_FINAL_SIZE)
kable(
  data.frame(
    Statistic = c(names(size_summary), "Std. dev."),
    `Hectares` = c(format(round(as.numeric(size_summary), 1), big.mark = ","),
                   format(round(sd(Fire_Disturbance_Point$FIRE_FINAL_SIZE), 1), big.mark = ",")),
    check.names = FALSE
  ),
  align = c("l", "r"), caption = "Final fire size (ha), all fires"
)

Final fire size (ha), all fires

Statistic	Hectares
Min.	0.0
1st Qu.	0.1
Median	0.2
Mean	147.3
3rd Qu.	1.0
Max.	191,810.6
Std. dev.	2,338.6

# Filter the data to fires between 15 and 99 hectares
fire_data_filtered <- fire_data %>%
  filter(FIRE_FINAL_SIZE > 15, FIRE_FINAL_SIZE <= 99)

ggplot(fire_data_filtered, aes(x = FIRE_FINAL_SIZE)) +
  geom_histogram(binwidth = 1, fill = ms_colours[["orange"]], colour = NA) +
  labs(title = "Distribution of final fire size",
       subtitle = "Fires between 15 and 99 hectares",
       x = "Fire final size (hectares)", y = "Number of fires") +
  theme_mikael()

Trends Over Time

Median fire size by year

median_fire_size_by_year <- fire_data %>%
  group_by(FIRE_YEAR) %>%
  filter(FIRE_FINAL_SIZE > 10, FIRE_FINAL_SIZE <= 100) %>%
  summarize(Median_Size = median(FIRE_FINAL_SIZE, na.rm = TRUE))

ggplot(median_fire_size_by_year, aes(x = FIRE_YEAR, y = Median_Size)) +
  geom_line(colour = ms_colours[["orange"]], linewidth = 0.8) +
  geom_point(colour = ms_colours[["orange"]], size = 1.8) +
  labs(title = "Median fire size by year",
       subtitle = "Fires between 10 and 100 hectares",
       x = "Year", y = "Median fire size (hectares)") +
  theme_mikael()

Average fire size by year, with linear trend

avg_fire_size_by_year <- fire_data %>%
  group_by(FIRE_YEAR) %>%
  filter(FIRE_FINAL_SIZE > 15, FIRE_FINAL_SIZE <= 100) %>%
  summarize(Average_Size = mean(FIRE_FINAL_SIZE, na.rm = TRUE))

# Linear regression on year
linear_model <- lm(Average_Size ~ FIRE_YEAR, data = avg_fire_size_by_year)
coef <- coef(linear_model)
summary_stats <- summary(linear_model)
r_squared <- summary_stats$r.squared
adj_r_squared <- summary_stats$adj.r.squared
f_value <- summary_stats$fstatistic[1]
p_value <- pf(f_value, df1 = summary_stats$fstatistic[2], df2 = summary_stats$fstatistic[3], lower.tail = FALSE)

stat_label <- paste(
  sprintf("y = %.2fx + %.2f", coef[2], coef[1]),
  sprintf("R² = %.3f", r_squared),
  sprintf("Adj. R² = %.3f", adj_r_squared),
  sprintf("F = %.2f", f_value),
  sprintf("p = %.3g", p_value),
  sep = "   ·   "
)

ggplot(avg_fire_size_by_year, aes(x = FIRE_YEAR, y = Average_Size)) +
  geom_line(colour = ms_colours[["orange"]], linewidth = 0.8) +
  geom_point(colour = ms_colours[["orange"]], size = 1.8) +
  geom_smooth(method = "lm", colour = ms_colours[["red"]], se = FALSE, linewidth = 0.8) +
  labs(title = "Average fire size by year, Ontario",
       subtitle = stat_label,
       x = "Year", y = "Average fire size (hectares)") +
  theme_mikael()

kable(
  data.frame(
    Term = c("Slope (ha per year)", "Intercept", "R²", "Adjusted R²", "F statistic", "p-value"),
    Value = c(sprintf("%.3f", coef[2]), sprintf("%.1f", coef[1]), sprintf("%.3f", r_squared),
              sprintf("%.3f", adj_r_squared), sprintf("%.2f", f_value), sprintf("%.3g", p_value))
  ),
  align = c("l", "r"), caption = "Linear model: average fire size (15–100 ha) ~ year"
)

Linear model: average fire size (15–100 ha) ~ year

Term	Value
Slope (ha per year)	0.150
Intercept	-258.4
R²	0.157
Adjusted R²	0.139
F statistic	8.41
p-value	0.00576

Fire frequency by year

fire_frequency_by_year <- fire_data %>%
  filter(FIRE_FINAL_SIZE > 15, FIRE_FINAL_SIZE <= 50) %>%
  group_by(FIRE_YEAR) %>%
  summarise(Frequency = n())

ggplot(fire_frequency_by_year, aes(x = FIRE_YEAR, y = Frequency)) +
  geom_col(fill = ms_colours[["orange"]], width = 0.75) +
  labs(title = "Fire frequency by year",
       subtitle = "Fires between 15 and 50 hectares",
       x = "Year", y = "Number of fires") +
  theme_mikael()

Where Fires Happen

ggplot() +
  base_map() +
  geom_sf(data = fires_sf, colour = ms_colours[["red"]], shape = ".") +
  labs(title = "Fire disturbance points, Ontario",
       subtitle = sprintf("All %s fires, %d–%d", format(nrow(fires_sf), big.mark = ","),
                          min(fires_sf$FIRE_YEAR), max(fires_sf$FIRE_YEAR))) +
  map_frame()

Causes

cause_counts <- Fire_Disturbance_Point %>%
  filter(FIRE_GENERAL_CAUSE_CODE %in% names(cause_labels)) %>%
  count(FIRE_GENERAL_CAUSE_CODE) %>%
  mutate(label = cause_labels[FIRE_GENERAL_CAUSE_CODE],
         label = reorder(label, n))

ggplot(cause_counts, aes(x = n, y = label, fill = FIRE_GENERAL_CAUSE_CODE)) +
  geom_col(width = 0.7, show.legend = FALSE) +
  geom_text(aes(label = format(n, big.mark = ",")), hjust = -0.15, colour = ms_text, size = 3.5) +
  scale_fill_manual(values = cause_colours) +
  scale_x_continuous(expand = expansion(mult = c(0, 0.12)), labels = scales::comma) +
  labs(title = "Fires by general cause", x = "Number of fires", y = NULL) +
  theme_mikael() +
  theme(panel.grid.major.y = element_blank())

ggplot() +
  base_map() +
  geom_sf(data = filter(fires_sf, FIRE_GENERAL_CAUSE_CODE %in% names(cause_labels)),
          aes(colour = FIRE_GENERAL_CAUSE_CODE), shape = ".") +
  scale_colour_manual(values = cause_colours, labels = cause_labels, name = "Cause") +
  guides(colour = guide_legend(override.aes = list(shape = 16, size = 3))) +
  labs(title = "Fire disturbance points by cause") +
  map_frame() +
  theme(legend.position = "right")

Large fires (over 40 hectares)

large_fires <- filter(fires_sf, FIRE_FINAL_SIZE > 40, FIRE_GENERAL_CAUSE_CODE %in% names(cause_labels))

ggplot() +
  base_map() +
  geom_sf(data = large_fires, aes(colour = FIRE_GENERAL_CAUSE_CODE), size = 0.5, alpha = 0.8) +
  scale_colour_manual(values = cause_colours, labels = cause_labels, name = "Cause") +
  guides(colour = guide_legend(override.aes = list(size = 3, alpha = 1))) +
  labs(title = "Fires over 40 hectares, by cause",
       subtitle = sprintf("%s fires", format(nrow(large_fires), big.mark = ","))) +
  map_frame() +
  theme(legend.position = "right")

large_lightning <- filter(fires_sf, FIRE_FINAL_SIZE > 40, FIRE_GENERAL_CAUSE_CODE == "LTG")

ggplot() +
  base_map() +
  geom_sf(data = large_lightning, colour = cause_colours[["LTG"]], size = 0.5, alpha = 0.85) +
  labs(title = "Lightning-caused fires over 40 hectares",
       subtitle = sprintf("%s fires", format(nrow(large_lightning), big.mark = ","))) +
  map_frame()

Interactive map of large fires

Fires over 40 hectares. Click a point for its year, cause and size.

large_ll <- st_transform(large_fires, 4326)
pal <- colorFactor(unname(cause_colours), levels = names(cause_colours))

leaflet(large_ll, width = "100%", height = 560) %>%
  addProviderTiles(providers$CartoDB.DarkMatter) %>%
  addCircleMarkers(
    radius = ~pmin(8, 2 + log10(FIRE_FINAL_SIZE)), stroke = FALSE, fillOpacity = 0.75,
    color = ~pal(FIRE_GENERAL_CAUSE_CODE),
    popup = ~sprintf("<b>%d</b><br>%s<br>%s ha", FIRE_YEAR, cause_labels[FIRE_GENERAL_CAUSE_CODE],
                     format(round(FIRE_FINAL_SIZE, 1), big.mark = ","))
  ) %>%
  addLegend("bottomright", colors = unname(cause_colours), labels = unname(cause_labels),
            title = "Cause", opacity = 1)

Clustering

K-means (k = 5), fires over 40 hectares

subset_data <- subset(Fire_Disturbance_Point, FIRE_FINAL_SIZE > 40)

k <- 5
set.seed(123)
subset_data$cluster <- kmeans(subset_data[, c("X", "Y")], centers = k)$cluster

ggplot(subset_data, aes(x = X, y = Y, colour = factor(cluster))) +
  geom_point(size = 0.5) +
  scale_colour_manual(values = unname(ms_colours[c("red", "orange", "green", "blue", "purple")]), name = "Cluster") +
  guides(colour = guide_legend(override.aes = list(size = 3))) +
  scale_x_continuous(breaks = seq(-95, -75, by = 2), limits = c(-95, -75)) +
  scale_y_continuous(breaks = seq(41, 57, by = 2), limits = c(41, 57)) +
  labs(title = "K-means clusters of fires over 40 hectares", x = "Longitude", y = "Latitude") +
  theme_mikael()

DBSCAN, lightning-caused fires over 40 hectares

DBSCAN (eps = 0.25°, minimum 10 points) finds dense groups of large lightning fires. The outlined boxes mark the main cluster regions.

eps_value <- 0.25
minPts_value <- 10
set.seed(123)
subset_data <- subset(Fire_Disturbance_Point, FIRE_FINAL_SIZE > 40 & FIRE_GENERAL_CAUSE_CODE == "LTG")
subset_data$cluster <- dbscan(subset_data[, c("X", "Y")], eps = eps_value, minPts = minPts_value)$cluster

# Cluster regions (lon/lat boxes)
boxes <- list(
  list(lon = c(-94.59020, -91.18305), lat = c(49.18557, 51.35367), col = "#d62221"),
  list(lon = c(-92.73902, -88.98409), lat = c(53.38905, 55.59270), col = "#8f3f71"),
  list(lon = c(-89.72745, -85.90076), lat = c(53.38143, 55.60380), col = "#b16286"),
  list(lon = c(-90.45155, -86.74749), lat = c(51.96440, 54.20702), col = "#98971a"),
  list(lon = c(-91.83242, -88.20501), lat = c(51.98362, 54.18996), col = "#d79921"),
  list(lon = c(-86.21592, -82.72439), lat = c(50.57883, 52.77495), col = "#458592")
)

# Noise points (cluster 0) in faint grey, clusters in distinct colours
n_clusters <- max(subset_data$cluster)
cluster_cols <- c("0" = "#8c857755", setNames(viridis(n_clusters, option = "H"), seq_len(n_clusters)))

p <- ggplot(subset_data, aes(x = X, y = Y, colour = factor(cluster))) +
  geom_point(size = 0.75) +
  scale_colour_manual(values = cluster_cols, guide = "none") +
  scale_x_continuous(breaks = seq(-95, -75, by = 2), limits = c(-95, -75)) +
  scale_y_continuous(breaks = seq(41, 57, by = 2), limits = c(41, 57)) +
  labs(title = "DBSCAN clusters, lightning fires over 40 hectares",
       subtitle = "Grey points are noise (not part of any cluster)", x = "Longitude", y = "Latitude") +
  theme_mikael()

for (b in boxes) {
  p <- p + annotate("rect", xmin = b$lon[1], xmax = b$lon[2], ymin = b$lat[1], ymax = b$lat[2],
                    fill = NA, colour = b$col, linewidth = 0.9)
}
p

All lightning-caused fires

ggplot() +
  base_map() +
  geom_sf(data = filter(fires_sf, FIRE_GENERAL_CAUSE_CODE == "LTG"),
          colour = cause_colours[["LTG"]], shape = ".") +
  labs(title = "All lightning-caused fires",
       subtitle = sprintf("%s fires", format(sum(fires_sf$FIRE_GENERAL_CAUSE_CODE == "LTG"), big.mark = ","))) +
  map_frame()

---

Data: Ontario Ministry of Natural Resources and Forestry, Fire Disturbance Point. Basemap: Natural Earth.