# ============================================================
# Knaappila A., Lindberg H., Tanskanen H., Aakala T. (2026). Ignition probability and fire size in Finnish boreal forests: a nationwide analysis of 10 years of rescue services data. Silva Fennica vol. 60 no. 3 article id 26006. 
# https://doi.org/10.14214/sf.26006

# Supplementary file S2. R-script.

# R. Core Team. R: A language and environment for statistical computing. Version 4.5.3.
# ============================================================

# Workflow:
# 1) Load required libraries
# 2) Read PRONTO and pseudo-absence data (spatial data already joined in ArcGIS Pro)
# 3) Rename PRONTO columns
# 4) Filter forest fires
# 5) Remove duplicate coordinates
# 6) Convert dates to day of year
# 7) Reclassify ignition causes
# 8) Remove unrealistic burnt area values
# 9) Calculate fire size in forested terrain types
# 10) Remove zeros and convert to hectares
# 11) Clean FFI
# 12) Fix canopy cover based on terrain type
# 13) Calculate dominant species
# 14) Label presence and pseudo-absence data
# 15) Remove rows with missing values and balance sample size
# 16) Convert variables to appropriate data types
# 17) Bind datasets
# 18) Fit ignition probability GAMs
# 19) Evaluate performance and effects (ignition probability)
# 20) Calculate spatial autocorrelation (ignition probability)
# 21) Fit fire size GAMs
# 22) Evaluate performance and effects (fire size)
# 23) Calculate spatial autocorrelation (fire size)


# ============================================================
# 1) Load required libraries
# ============================================================
library(dplyr)
library(mgcv)
library(sf)
library(spdep)
library(stringr)
library(ggplot2)
library(visreg)

# ============================================================
# 2) Read PRONTO and pseudo-absence data (spatial data already joined in ArcGIS Pro)
# ============================================================
data <- read.csv("path_to_pronto_data.csv")
ctrl <- read.csv("path_to_absence_data.csv")


# ============================================================
# 3) Rename PRONTO columns
# ============================================================
data$year <- data$Vuosi
data$date <- data$Ilmoitusaika
data$easting <- data$`Itäkoord. (7 num.)`
data$northing <- data$`Pohjoiskoord. (7 num.)`
data$accident_type <- data$`Onnettomuustyyppi 1`
data$terrain_type_1 <- data$`Maastotyyppi 1`
data$terrain_type_2 <- data$`Maastotyyppi 2`
data$ffi <- data$Metsäpaloindeksi
data$burnt_arrival <- data$`Palanut maastoala palokunnan saapuessa (m2)`
data$burnt1 <- data$`Palanut maastoala tilanteen lopussa maastotyypissä 1 (m2)`
data$burnt2 <- data$`Palanut maastoala tilanteen lopussa maastotyypissä 2 (m2)`
data$initiator <- data$`Arvio, mikä aiheutti tulipalon`
data$ignition_cause <- data$`Arvio tulipalon syttymissyystä (lisätty 1.1.2009)`


# ============================================================
# 4) Filter forest fires
# ============================================================
data <- data[data$accident_type == "Maastopalo", ]
data <- data[data$terrain_type_1 == "Hakkuualue, raivio tai metsäaukea" 
             | data$terrain_type_1 == "Vars. metsä, talousmetsä, aluskasvillisuus metsäalueella"
             | data$terrain_type_1 == "Muu suoalue", ]


# ============================================================
# 5) Remove duplicate coordinates
# ============================================================
data <- (data %>%
           group_by(easting, northing) %>%
           filter(n() == 1))


# ============================================================
# 6) Convert dates to day of year
# ============================================================
data$doy <- strftime(data$date, format = "%j")


# ============================================================
# 7) Reclassify ignition causes
# ============================================================
data$cause1 <- dplyr::recode(data$initiator,
                             "Ihmisen toiminta"="1",
                             "Ei voida arvioida"="3",
                             "Luonnontapahtuma tai -ilmiö"="2",
                             "Koneen tai laitteen vika"="1",
                             "Muu syy"="3",
                             "Palovaarallinen aine"="1",      
                             "Eläin"="3",
                             .default = "3",
                             .missing = "3")

data$cause2 <- dplyr::recode(data$ignition_cause,
                             "Nuotio, grilli"="1", 
                             "Kulotus"="1",                                                           
                             "Tulitikku, muu tulentekoväline"="1", 
                             "Ei voida arvioida"="3",
                             "Salama"="2", 
                             "Lasten tulen käsittely"="1",
                             "Muu syy"="3", 
                             "Roskien, risujen poltto"="1",
                             "Uudelleensyttyminen"="3", 
                             "Mekaaninen kipinä, iskukipinä"="1",
                             "Muu luonnonsyy"="3", 
                             "Kuuma tai hehkuva esine tai tuhka"="1",
                             "Savuke tai muu tupakka-aine"="1", 
                             "Tahallaan sytytetty palo"="1",
                             "Kipinä tai kekäle tulisijasta tai hormista"="1", 
                             "Ilotulite, pyrotekniset tuotteet"="1",
                             "Koneen tai laitteen väärä käyttö"="1", 
                             "Räjähdys"="3",
                             "Kipinä koneesta tai laitteesta"="1", 
                             "Koneen tai laitteen vika, häiriö tai huollon laiminlyönti"="1",
                             "Muu ruoanvalmistus"="1", 
                             "Sähkölaitteen tai -asennuksen vika, häiriö tai huollon laiminlyönti"="1",
                             "Itsesyttymä"="3", 
                             "Hankauslämpö"="1",
                             "Kynttilä, tuikku, soihtu, roihu"="1", 
                             "Tulityö"="1",
                             "Riittämätön suojaetäisyys"="1", 
                             "Valvomaton ruoanvalmistus"="1",
                             "Vaurio tulisijassa tai hormissa"="1", 
                             "Tuotantoprosessin häiriö"="1",
                             "Nuotio"="1", 
                             "Grilli"="1",
                             .default = "3",
                             .missing = "3")

data$ignition <- dplyr::case_when(
  data$cause1 %in% c ("1") & 
    data$cause2 %in% c ("1") ~ "Human",
  data$cause1 %in% c ("1") & 
    data$cause2 %in% c ("2") ~ "Unknown",
  data$cause1 %in% c ("1") & 
    data$cause2 %in% c ("3") ~ "Human",
  data$cause1 %in% c ("2") & 
    data$cause2 %in% c ("1") ~ "Human",
  data$cause1 %in% c ("2") & 
    data$cause2 %in% c ("2") ~ "Lightning",
  data$cause1 %in% c ("2") & 
    data$cause2 %in% c ("3") ~ "Unknown",
  data$cause1 %in% c ("3") & 
    data$cause2 %in% c ("1") ~ "Human",
  data$cause1 %in% c ("3") &
    data$cause2 %in% c ("2") ~ "Unknown",
  data$cause1 %in% c ("3") & 
    data$cause2 %in% c ("3") ~ "Unknown")


# ============================================================
# 8) Remove unrealistic burnt area values
# ============================================================
data$burnt1[is.na(data$burnt1)] <- 0
data$burnt2[is.na(data$burnt2)] <- 0
data$burnt_arrival[is.na(data$burnt_arrival)] <- 0

data$burnt_total <- with(data, data$burnt1 + data$burnt2)
data$burnt_ratio <- with(data, burnt_total / burnt_arrival)

data <- data[data$burnt_total!= 0, ]
data <- data[data$burnt_ratio != 10000, ]
data <- data[data$burnt_ratio != 1000, ]
data <- data[data$burnt_ratio != 100, ]
data <- data[data$burnt_ratio != 10, ]
data <- data[data$burnt_ratio != 0, ]
data <- data[!(data$burnt_ratio < 1.0), ]


# ============================================================
# 9) Calculate fire size in forested terrain types
# ============================================================
data$burnt_forest1 <- ifelse(data$terrain_type_1 %in% "Vars. metsä, talousmetsä, aluskasvillisuus metsäalueella" |
                               data$terrain_type_1 %in% "Hakkuualue, raivio tai metsäaukea" |
                               data$terrain_type_1 %in% "Muu suoalue",
                             data$burnt1, 0)
data$burnt_forest2 <- ifelse(data$ terrain_type_2 %in% "Vars. metsä, talousmetsä, aluskasvillisuus metsäalueella" |
                               data$ terrain_type_2 %in% "Hakkuualue, raivio tai metsäaukea" |
                               data$ terrain_type_2 %in% "Muu suoalue",
                             data$burnt2, 0)
data$fire_size <- data$burnt_forest1 + data$burnt_forest2


# ============================================================
# 10) Remove zeros and convert to hectares
# ============================================================
data <- data[data$fire_size!= 0, ]
data$fire_size <- data$fire_size/10000


# ============================================================
# 11) Clean FFI
# ============================================================
data <- data %>%
  mutate(
    ffi = str_trim(ffi),
    ffi = str_replace_all(ffi, ",", "."),
    ffi = na_if(ffi, ""),
    ffi = as.numeric(ffi),
    ffi = ifelse(ffi >= 1 & ffi <= 6,
                 round(ffi, 1),
                 NA))


# ============================================================
# 12) Fix canopy cover based on terrain type
# ============================================================
data$cover <- ifelse(data$terrain_type_1 == "Hakkuualue, raivio tai metsäaukea",
                     0, data$cover)


# ============================================================
# 13) Calculate dominant species
# ============================================================
data$volume_total <- with(data, pine + spruce + birch + otherdec)
data$volume_deciduous <- with(data, birch + otherdec)
data$share_deciduous <- with(data, volume_deciduous / volume_total)
data$share_pine <- with(data, pine / volume_total)
data$share_spruce <- with(data, spruce / volume_total)

data$species <- ifelse(data$share_deciduous > 0.5, "Deciduous",
                       ifelse(data$share_pine > 0.5, "Pine",
                              ifelse(data$share_spruce > 0.5, "Spruce", "Mixed")))

ctrl$volume_total <- with(ctrl, pine + spruce + birch + otherdec)
ctrl$volume_deciduous <- with(ctrl, birch + otherdec)
ctrl$share_deciduous <- with(ctrl, volume_deciduous / volume_total)
ctrl$share_pine <- with(ctrl, pine / volume_total)
ctrl$share_spruce <- with(ctrl, spruce / volume_total)

ctrl$species <- ifelse(ctrl$share_deciduous > 0.5, "Deciduous",
                       ifelse(ctrl$share_pine > 0.5, "Pine",
                              ifelse(ctrl$share_spruce > 0.5, "Spruce", "Mixed")))


# ============================================================
# 14) Label presence and pseudo-absence data
# ============================================================
data$fire <- 1
ctrl$fire <- 0

data <- data %>%
  dplyr::select(c(year, easting, northing, fire, fire_size,
                  trafficability, slope, cover, species, doy, ffi, ignition, population, road))

ctrl <- ctrl %>%
  dplyr::select(year, easting, northing, fire,
                trafficability, slope, cover, species, doy, population, road)

# ============================================================
# 15) Remove rows with missing values and balance sample size
# ============================================================
data <- data[complete.cases(data),]
ctrl <- ctrl[complete.cases(ctrl),]

set.seed(123)
ctrl <- ctrl[sample(nrow(ctrl), nrow(data)), ]


# ============================================================
# 16) Convert variables to appropriate data types
# ============================================================
data$year <- as.integer(data$year)
data$easting <- as.integer(data$easting)
data$northing <- as.integer(data$northing)
data$fire <- as.numeric(data$fire)
data$fire_size <- as.numeric(data$fire_size)
data$trafficability <- as.factor(data$trafficability)
data$slope <- as.numeric(data$slope)
data$cover <- as.numeric(data$cover)
data$species <- as.factor(data$species)
data$doy <- as.integer(data$doy)
data$ffi <- as.numeric(data$ffi)
data$ignition <- as.factor(data$ignition)
data$population <- as.numeric(data$population)
data$road <- as.numeric(data$road)

ctrl$year <- as.integer(ctrl$year)
ctrl$easting <- as.integer(ctrl$easting)
ctrl$northing <- as.integer(ctrl$northing)
ctrl$fire <- as.numeric(ctrl$fire)
ctrl$trafficability <- as.factor(ctrl$trafficability)
ctrl$slope <- as.numeric(ctrl$slope)
ctrl$cover <- as.numeric(ctrl$cover)
ctrl$species <- as.factor(ctrl$species)
ctrl$doy <- as.integer(ctrl$doy)
ctrl$population <- as.numeric(ctrl$population)
ctrl$road <- as.numeric(ctrl$road)


# ============================================================
# 17) Bind datasets
# ============================================================
data1 <- bind_rows(data, ctrl)
data2 <- data1[data1$fire == 1,]


# ============================================================
# 18) Fit ignition probability GAMs
# ============================================================
gam_ignition_1 <- gam(fire ~
                        trafficability +
                        s(slope, k=3) +
                        s(cover) +
                        species +
                        s(doy) +
                        s(population) +
                        road,
                      family = binomial(),
                      select = TRUE,
                      data = data1)

gam_ignition_2 <- gam(fire ~
                        trafficability +
                        s(slope, k=3) +
                        s(cover) +
                        species +
                        s(doy) +
                        s(population) +
                        road +
                        s(easting, northing),
                      family = binomial(),
                      select = TRUE,
                      data = data1)


# ============================================================
# 19) Evaluate performance and effects (ignition probability)
# ============================================================
summary(gam_ignition_1)
anova(gam_ignition_1)
visreg(gam_ignition_1, rug = TRUE, scale = "response")
round(concurvity(gam_ignition_1, full = TRUE), 2)
mean((predict(gam_ignition_1, type = "response") >= 0.5) == gam_ignition_1$model[[1]]) * 100

summary(gam_ignition_2)
anova(gam_ignition_2)
visreg(gam_ignition_2, rug = TRUE, scale = "response")
round(concurvity(gam_ignition_2, full = TRUE), 2)
mean((predict(gam_ignition_2, type = "response") >= 0.5) == gam_ignition_2$model[[1]]) * 100

# ============================================================
# 20) Calculate spatial autocorrelation (ignition probability)
# ============================================================
data1$resid1 <- residuals(gam_ignition_1, type = "pearson")
data1$resid2 <- residuals(gam_ignition_2, type = "pearson")
data1_spatial <- st_as_sf(data1, coords = c("easting", "northing"), crs = 2393)
distances1 <- seq(1000, 50000, by = 5000)
coords1 <- st_coordinates(data1_spatial)

moran_ignition_1 <- do.call(rbind, lapply(distances1, function(d) {
  cat("Distance =", d/1000, "km\n")
  flush.console()
  nb <- dnearneigh(coords1, d1 = 0, d2 = d)
  lw <- nb2listw(nb, zero.policy = TRUE)
  mt <- moran.test(
    data1_spatial$resid1,
    lw,
    zero.policy = TRUE)
  data.frame(
    Distance = d,
    MoransI = unname(mt$estimate["Moran I statistic"]),
    Pvalue = mt$p.value)
}))

moran_ignition_2 <- do.call(rbind, lapply(distances1, function(d) {
  cat("Distance =", d/1000, "km\n")
  flush.console()
  nb <- dnearneigh(coords1, d1 = 0, d2 = d)
  lw <- nb2listw(nb, zero.policy = TRUE)
  mt <- moran.test(
    data1_spatial$resid2,
    lw,
    zero.policy = TRUE)
  data.frame(
    Distance = d,
    MoransI = unname(mt$estimate["Moran I statistic"]),
    Pvalue = mt$p.value)
}))

moran_ignition_1$model  <- "Non-spatial GAM"
moran_ignition_2$model <- "Spatial GAM"
moran_ignition <- rbind(moran_ignition_1, moran_ignition_2)

ggplot(moran_ignition,
       aes(x = Distance / 1000,
           y = MoransI,
           linetype = model)) +
  geom_hline(yintercept = 0) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  labs(x = "Distance (km)",
       y = "Moran's I",
       linetype = "") +
  theme_bw()



# ============================================================
# 21) Fit fire size GAMs
# ============================================================
gam_size_1 <- gam(log(fire_size) ~
                    s(doy) +
                    s(ffi) +
                    ignition +
                    slope +
                    trafficability +
                    s(cover) +
                    species +
                    s(population) +
                    road,
                  family = gaussian(),
                  select = TRUE,
                  data = data2)

gam_size_2 <- gam(log(fire_size) ~
                    s(doy) +
                    s(ffi) +
                    ignition +
                    slope +
                    trafficability +
                    s(cover) +
                    species +
                    s(population) +
                    road +
                    s(easting, northing),
                  family = gaussian(),
                  select = TRUE,
                  data = data2)


# ============================================================
# 22) Evaluate performance and effects (fire size)
# ============================================================
summary(gam_size_1)
anova(gam_size_1)
visreg(gam_size_1, rug = TRUE, trans = exp)
round(concurvity(gam_size_1, full = TRUE), 2)

summary(gam_size_2)
anova(gam_size_2)
visreg(gam_size_2, rug = TRUE, trans = exp)
round(concurvity(gam_size_2, full = TRUE), 2)


# ============================================================
# 23) Calculate spatial autocorrelation (fire size)
# ============================================================
data2$resid1 <- residuals(gam_size_1, type = "pearson")
data2$resid2 <- residuals(gam_size_2, type = "pearson")
data2_spatial <- st_as_sf(data2, coords = c("easting", "northing"), crs = 2393)
distances2 <- seq(1000, 50000, by = 5000)
coords2 <- st_coordinates(data2_spatial)

moran_size_1 <- do.call(rbind, lapply(distances2, function(d) {
  cat("Distance =", d/1000, "km\n")
  flush.console()
  nb <- dnearneigh(coords2, d1 = 0, d2 = d)
  lw <- nb2listw(nb, zero.policy = TRUE)
  mt <- moran.test(
    data2_spatial$resid1,
    lw,
    zero.policy = TRUE)
  data.frame(
    Distance = d,
    MoransI = unname(mt$estimate["Moran I statistic"]),
    Pvalue = mt$p.value)
}))

moran_size_2 <- do.call(rbind, lapply(distances2, function(d) {
  cat("Distance =", d/1000, "km\n")
  flush.console()
  nb <- dnearneigh(coords2, d1 = 0, d2 = d)
  lw <- nb2listw(nb, zero.policy = TRUE)
  mt <- moran.test(
    data2_spatial$resid2,
    lw,
    zero.policy = TRUE)
  data.frame(
    Distance = d,
    MoransI = unname(mt$estimate["Moran I statistic"]),
    Pvalue = mt$p.value)
}))

moran_size_1$model  <- "Non-spatial GAM"
moran_size_2$model <- "Spatial GAM"
moran_size <- rbind(moran_size_1, moran_size_2)

ggplot(moran_size,
       aes(x = Distance / 1000,
           y = MoransI,
           linetype = model)) +
  geom_hline(yintercept = 0) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  labs(x = "Distance (km)",
       y = "Moran's I",
       linetype = "") +
  theme_bw()