## Reading layer `fire_points' from data source 
##   `https://drive.google.com/uc?export=download&id=16m9HJAn5uzSyiUdhdhoboPRRi9dDDQaT' 
##   using driver `GeoJSON'
## Simple feature collection with 1138 features and 33 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 2529244 ymin: 345248.9 xmax: 2680530 ymax: 484827
## Projected CRS: NAD83 / North Carolina (ftUS)
###THIS TURNS THE DATA INTO A POINT PATTERN OBJECT WITH A DEFINED BOUNDARY###
### USE THIS FOR PPA####
combined_wildfire <- st_union(wildfire_points)
bb_wildfire <- st_convex_hull(combined_wildfire) |> as.owin()
wildfire.ppp <- as.ppp(st_coordinates(wildfire_points), W = bb_wildfire)
# ---- Map with basemap ----
library(tmap)
tmap_mode("view")

tm_shape(wildfire_points) +
  tm_dots(col = "red", size = 0.05, alpha = 0.7)
# ---- Quadrat analysis (global structure) ----
q_count_wildfire <- quadratcount(wildfire.ppp, nx = 8, ny = 8)

plot(wildfire.ppp, main = "Wildfire Points with Quadrat Counts", pch = 20)
plot(q_count_wildfire, add = TRUE, col = "red")

# variance-to-mean ratio
var(as.vector(q_count_wildfire)) / mean(as.vector(q_count_wildfire))
## [1] 29.42872
# Monte Carlo quadrat test
q_test_wildfire <- quadrat.test(wildfire.ppp, nx = 8, ny = 8,
                                 method = "MonteCarlo", nsim = 99)
q_test_wildfire
## 
##  Conditional Monte Carlo test of CSR using quadrat counts
##  Test statistic: Pearson X2 statistic
## 
## data:  wildfire.ppp
## X2 = 1083.4, p-value = 0.02
## alternative hypothesis: two.sided
## 
## Quadrats: 54 tiles (irregular windows)
plot(q_test_wildfire)

# ---- Kernel Density Estimation ----
kernel_density_wildfire <- density.ppp(wildfire.ppp)
plot(kernel_density_wildfire, main = "KDE of Wildfire Points")

# ---- ANN with Monte Carlo simulations ----
# NOTE: Look at your quadrat.test() p-value above first.
# If p < 0.05 (reject CSR / evidence of a global trend) -> use the
#   INHOMOGENEOUS block (simulates from the KDE surface).
# If p >= 0.05 (fail to reject CSR / no clear global trend) -> use the
#   HOMOGENEOUS block (simulates CSR at a constant intensity).
# Only run ONE of the two blocks below - delete/comment out the other.

## -- HOMOGENEOUS version --
ANN_obs_wildfire <- mean(nndist(wildfire.ppp))
nsim <- 99
ANN_sim_wildfire <- numeric(nsim)
lambda_wildfire <- wildfire.ppp$n / area.owin(wildfire.ppp$window)

for (i in 1:nsim) {
  sim <- rpoispp(lambda_wildfire, win = wildfire.ppp$window)
  ANN_sim_wildfire[i] <- mean(nndist(sim))
}

N_greater <- sum(ANN_sim_wildfire >= ANN_obs_wildfire)
p_value_ann <- 2 * min((N_greater + 1) / (nsim + 1),
                        (nsim - N_greater + 1) / (nsim + 1))
p_value_ann
## [1] 0.02
hist(ANN_sim_wildfire,
     main = "Monte Carlo Distribution of ANN (Wildfire, Homogeneous)",
     xlab = "Average Nearest Neighbor Distance")
abline(v = ANN_obs_wildfire, col = "red", lwd = 2)

## -- INHOMOGENEOUS version --
ANN_obs_wildfire <- mean(nndist(wildfire.ppp))
nsim <- 99
ANN_sim_wildfire <- numeric(nsim)

for (i in 1:nsim) {
  sim <- rpoispp(kernel_density_wildfire)
  ANN_sim_wildfire[i] <- mean(nndist(sim))
}

N_greater <- sum(ANN_sim_wildfire >= ANN_obs_wildfire)
p_value_ann <- 2 * min((N_greater + 1) / (nsim + 1),
                        (nsim - N_greater + 1) / (nsim + 1))
p_value_ann
## [1] 0.02
hist(ANN_sim_wildfire,
     main = "Monte Carlo Distribution of ANN (Wildfire, Inhomogeneous)",
     xlab = "Average Nearest Neighbor Distance")
abline(v = ANN_obs_wildfire, col = "red", lwd = 2)

# ---- L-function and envelope (border correction, not iso) ----
# Same logic as above: choose Lest (homogeneous) or Linhom (inhomogeneous)
# based on your quadrat test result. Only run ONE block.

## -- HOMOGENEOUS version --
l_funct_wildfire <- envelope(
  wildfire.ppp,
  fun = Lest,
  nsim = 99,
  funargs = list(correction = "border"),
  transform = expression(. - r)
)
## Generating 99 simulated realisations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
## 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
## 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 
## 99.
## 
## Done.
plot(l_funct_wildfire, legend = FALSE)

## -- INHOMOGENEOUS version --
l_funct_wildfire <- envelope(
  wildfire.ppp,
  fun = Linhom,
  nsim = 99,
  funargs = list(correction = "border"),
  transform = expression(. - r)
)
## Generating 99 simulated realisations of CSR  ...
## 1,
## 2,
## 3,
## 4,
## 5,
## 6,
## 7,
## 8,
## 9,
## 10,
## 11,
## 12,
## 13,
## 14,
## 15,
## 16,
## 17,
## 18,
## 19,
## 20,
## 21,
## 22,
## 23,
## 24,
## 25,
## 26,
## 27,
## 28,
## 29,
## 30,
## 31,
## 32,
## 33,
## 34,
## 35,
## 36,
## 37,
## 38,
## 39,
## 40,
## 41,
## 42,
## 43,
## 44,
## 45,
## 46,
## 47,
## 48,
## 49,
## 50,
## 51,
## 52,
## 53,
## 54,
## 55,
## 56,
## 57,
## 58,
## 59,
## 60,
## 61,
## 62,
## 63,
## 64,
## 65,
## 66,
## 67,
## 68,
## 69,
## 70,
## 71,
## 72,
## 73,
## 74,
## 75,
## 76,
## 77,
## 78,
## 79,
## 80,
## 81,
## 82,
## 83,
## 84,
## 85,
## 86,
## 87,
## 88,
## 89,
## 90,
## 91,
## 92,
## 93,
## 94,
## 95,
## 96,
## 97,
## 98,
## 
## 99.
## 
## Done.
plot(l_funct_wildfire, legend = FALSE)

Interpretation The evidence points to both first-order and second-order processes shaping this wildfire pattern, not just one or the other.

First-order evidence comes from the KDE surface, which shows a clear non-random gradient in fire density across the study area rather than a flat, uniform intensity. This is the signature of an underlying environmental or land-use trend influencing where fires are more or less likely to occur overall. The very high variance-to-mean ratio (29.43) from the quadrat analysis reinforces this: counts vary far more across space than a homogeneous Poisson process would allow.

Second-order evidence comes from the fact that clustering persists even after that global trend is accounted for. The ANN test simulated against the inhomogeneous (KDE-based) intensity surface still returned a highly significant result (p = 0.02, with every one of the 99 simulations showing a larger nearest-neighbor distance than the observed data). This means fires sit closer to each other than the underlying density surface alone would predict. Likewise, the inhomogeneous L-function remains above the simulation envelope across most distances, confirming points are attracting/co-occurring with each other locally, beyond what the broader spatial trend explains.

Taken together, this is a combination pattern: a first-order gradient (driving where fire risk is generally higher) layered with genuine second-order clustering (individual fires occurring near one another beyond what that gradient predicts).

Hypothesis: In North Carolina, a large share of wildfires are historically linked to escaped human-caused ignitions. Particularly debris and agricultural burning that are concentrated near roads and rural development at the wildland–urban interface. This could plausibly explain both patterns observed here: the first-order gradient may reflect proximity to roads/development (more potential ignition sources in accessible areas), while the second-order clustering could reflect multiple escaped burns occurring in the same neighborhoods or along the same rural corridors during shared dry, windy conditions favorable to fire spread.