Tutorial

Point pattern analysis helps separate two kinds of spatial processes. A first-order process reflects variation in the surrounding environment, while a second-order process reflects attraction or repulsion among events themselves.

Reading in the simulated data

data_candidates <- c(
  "point_patterns.RData",
  file.path("upload", "point_patterns.RData")
)

data_path <- data_candidates[file.exists(data_candidates)][1]

if (is.na(data_path)) {
  stop("Place point_patterns.RData in the same folder as this Rmd file.")
}

load(data_path)

Global structure: quadrat analysis

Sample pattern 1

plot(sample_pattern_1, main = "sample_pattern_1")

q_count1 <- quadratcount(sample_pattern_1, nx = 3, ny = 3)
plot(q_count1, main = "Quadrat counts: sample_pattern_1")

vmr1 <- var(as.vector(q_count1)) / mean(as.vector(q_count1))
vmr1
## [1] 5.595339
q_test1 <- quadrat.test(
  sample_pattern_1,
  nx = 3,
  ny = 3,
  method = "MonteCarlo",
  nsim = 99
)
q_test1
## 
##  Conditional Monte Carlo test of CSR using quadrat counts
##  Test statistic: Pearson X2 statistic
## 
## data:  sample_pattern_1
## X2 = 44.763, p-value = 0.02
## alternative hypothesis: two.sided
## 
## Quadrats: 3 by 3 grid of tiles
plot(q_test1, main = "Quadrat test: sample_pattern_1")

Q2 response. The bottom-right quadrat has the largest positive residual, while the top-left quadrat has the largest negative residual. In other words, there are many more points than expected on the right side of the study area and far fewer than expected on the upper-left side. This uneven east-west pattern is evidence of clustering caused by a broad first-order trend rather than complete spatial randomness.

The variance-to-mean ratio is 5.6. Because it is well above 1 and the Monte Carlo quadrat test is significant, sample pattern 1 is globally clustered.

Sample pattern 2

plot(sample_pattern_2, main = "sample_pattern_2")

q_count2 <- quadratcount(sample_pattern_2, nx = 3, ny = 3)
plot(q_count2, main = "Quadrat counts: sample_pattern_2")

vmr2 <- var(as.vector(q_count2)) / mean(as.vector(q_count2))
vmr2
## [1] 1.868421
q_test2 <- quadrat.test(
  sample_pattern_2,
  nx = 3,
  ny = 3,
  method = "MonteCarlo",
  nsim = 999
)
q_test2
## 
##  Conditional Monte Carlo test of CSR using quadrat counts
##  Test statistic: Pearson X2 statistic
## 
## data:  sample_pattern_2
## X2 = 14.947, p-value = 0.122
## alternative hypothesis: two.sided
## 
## Quadrats: 3 by 3 grid of tiles
plot(q_test2, main = "Quadrat test: sample_pattern_2")

Sample pattern 2 has a smaller variance-to-mean ratio (1.87) and its quadrat test is not significant. I therefore fail to reject complete spatial randomness at the global scale.

Global structure: kernel density

kernel_density1 <- density.ppp(sample_pattern_1)
plot(kernel_density1, main = "Kernel density: sample_pattern_1")

kernel_density2 <- density.ppp(sample_pattern_2)
plot(kernel_density2, main = "Kernel density: sample_pattern_2")

The first KDE shows a broad increase in intensity toward the right side. The second has several local peaks but does not show the same statistically significant global trend in the quadrat test.

Local structure: average nearest neighbor

nsim <- 999

ANN_obs1 <- mean(nndist(sample_pattern_1))
ANN_sim1 <- numeric(nsim)

for (i in seq_len(nsim)) {
  sim <- rpoispp(kernel_density1)
  ANN_sim1[i] <- mean(nndist(sim))
}

N_greater1 <- sum(ANN_sim1 >= ANN_obs1)
p_value1 <- 2 * min(
  (N_greater1 + 1) / (nsim + 1),
  (nsim - N_greater1 + 1) / (nsim + 1)
)
p_value1
## [1] 0.966
hist(
  ANN_sim1,
  main = "Monte Carlo ANN: sample_pattern_1",
  xlab = "Average nearest-neighbor distance"
)
abline(v = ANN_obs1, col = "red", lwd = 2)

ANN_obs2 <- mean(nndist(sample_pattern_2))
ANN_sim2 <- numeric(nsim)
lambda2 <- sample_pattern_2$n / area.owin(sample_pattern_2$window)

for (i in seq_len(nsim)) {
  sim <- rpoispp(lambda2, win = sample_pattern_2$window)
  ANN_sim2[i] <- mean(nndist(sim))
}

N_greater2 <- sum(ANN_sim2 >= ANN_obs2)
p_value2 <- 2 * min(
  (N_greater2 + 1) / (nsim + 1),
  (nsim - N_greater2 + 1) / (nsim + 1)
)
p_value2
## [1] 0.896
hist(
  ANN_sim2,
  main = "Monte Carlo ANN: sample_pattern_2",
  xlab = "Average nearest-neighbor distance"
)
abline(v = ANN_obs2, col = "red", lwd = 2)

Neither ANN test is significant, so these data do not provide evidence of very local clustering or dispersion beyond their selected Poisson reference processes.

Local structure: L-function

l_funct_pattern1 <- envelope(
  sample_pattern_1,
  fun = Linhom,
  nsim = 99,
  funargs = list(correction = "iso"),
  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_pattern1, legend = FALSE, main = "Inhomogeneous L: pattern 1")

l_funct_pattern2 <- envelope(
  sample_pattern_2,
  fun = Lest,
  nsim = 99,
  funargs = list(correction = "iso"),
  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_pattern2, legend = FALSE, main = "Homogeneous L: pattern 2")

For pattern 1, the inhomogeneous L-function stays near the simulation envelope after the global intensity trend is accounted for. For pattern 2, the observed line rises above the homogeneous envelope at intermediate distances, showing local clustering at those scales.

Mini-Challenge: Wildfires

Read and prepare the data

fire_url <- paste0(
  "https://drive.google.com/uc?export=download&id=",
  "16m9HJAn5uzSyiUdhdhoboPRRi9dDDQaT"
)

fire_candidates <- c("fire_points.geojson", "wildfire_points.geojson")
fire_local <- fire_candidates[file.exists(fire_candidates)][1]
fire_source <- if (is.na(fire_local)) fire_url else fire_local

wildfire_points <- st_read(fire_source, quiet = TRUE) |>
  st_transform(crs = 2264)

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)

wildfire.ppp
## Planar point pattern: 1138 points
## window: polygonal boundary
## enclosing rectangle: [2529243.8, 2680530.2] x [345248.9, 484827] units

The data include 1138 wildfire records in and around the Croatan National Forest area of coastal North Carolina.

Map with basemap

tmap_mode("view")

wildfire_wgs84 <- st_transform(wildfire_points, 4326)

# tmap changed its aesthetic names in version 4. This keeps the map
# compatible with both the older and newer syntax.
if (packageVersion("tmap") >= "4.0.0") {
  wildfire_map <- tm_shape(wildfire_wgs84) +
    tm_basemap("OpenStreetMap") +
    tm_dots(fill = "firebrick", size = 0.025, fill_alpha = 0.55)
} else {
  wildfire_map <- tm_shape(wildfire_wgs84) +
    tm_basemap("OpenStreetMap") +
    tm_dots(col = "firebrick", size = 0.025, alpha = 0.55)
}

wildfire_map
tmap_mode("plot")

The points are spread across the study area, but several visible concentrations occur in the center, south-central portion, and eastern side.

Quadrat analysis

set.seed(123)

q_count_fire <- quadratcount(wildfire.ppp, nx = 8, ny = 8)
plot(q_count_fire, main = "8 x 8 wildfire quadrat counts")

vmr_fire <- var(as.vector(q_count_fire)) / mean(as.vector(q_count_fire))
vmr_fire
## [1] 29.42872
q_test_fire <- quadrat.test(
  wildfire.ppp,
  nx = 8,
  ny = 8,
  method = "MonteCarlo",
  nsim = 99
)
q_test_fire
## 
##  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_fire, main = "Wildfire quadrat test")

The variance-to-mean ratio is 29.43, far above the CSR expectation of 1. The Monte Carlo p-value is 0.02, so I reject complete spatial randomness. The strong differences among quadrats show a first-order change in wildfire intensity across the study area.

Kernel density

# Scott's rule is stable even though this historical dataset contains
# repeated coordinates for fires recorded at the same location.
fire_bw <- bw.scott.iso(wildfire.ppp)
kernel_density_fire <- density.ppp(
  wildfire.ppp,
  sigma = fire_bw,
  edge = TRUE
)

fire_bw
##    sigma 
## 9855.057
plot(kernel_density_fire, main = "Kernel density of wildfire locations")
plot(wildfire.ppp, add = TRUE, pch = 16, cex = 0.25, cols = "white")

The KDE confirms that wildfire intensity is not constant. The largest hotspot is in the south-central part of the study area, with another broad high-intensity band extending northward and eastward.

Average nearest neighbor analysis

Because the quadrat analysis showed global inhomogeneity, the ANN simulations use the KDE intensity surface rather than a homogeneous CSR process.

set.seed(123)

nsim_fire <- 99
ANN_obs_fire <- mean(nndist(wildfire.ppp))
ANN_sim_fire <- numeric(nsim_fire)

for (i in seq_len(nsim_fire)) {
  sim <- rpoispp(kernel_density_fire)
  ANN_sim_fire[i] <- mean(nndist(sim))
}

N_greater_fire <- sum(ANN_sim_fire >= ANN_obs_fire)
p_value_fire <- 2 * min(
  (N_greater_fire + 1) / (nsim_fire + 1),
  (nsim_fire - N_greater_fire + 1) / (nsim_fire + 1)
)

ANN_obs_fire
## [1] 767.1578
mean(ANN_sim_fire)
## [1] 1655.483
p_value_fire
## [1] 0.02
hist(
  ANN_sim_fire,
  main = "Monte Carlo ANN for wildfire locations",
  xlab = "Average nearest-neighbor distance (feet)"
)
abline(v = ANN_obs_fire, col = "red", lwd = 2)

The observed ANN is much smaller than the simulated values and the empirical two-sided p-value is 0.02. This is evidence of very local clustering beyond the broad intensity surface.

Inhomogeneous L-function

The global test was significant, so I use Linhom and simulate from the KDE intensity. The convex-hull window is irregular, so the edge correction is "border" as instructed.

set.seed(123)

l_funct_fire <- envelope(
  wildfire.ppp,
  fun = Linhom,
  simulate = expression(rpoispp(kernel_density_fire)),
  nsim = 99,
  funargs = list(
    lambda = kernel_density_fire,
    correction = "border"
  ),
  transform = expression(. - r)
)
## Generating 99 simulations by evaluating expression  ...
## 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_fire,
  legend = FALSE,
  main = "Inhomogeneous L-function: wildfire locations"
)

The observed function rises above the simulation envelope at short distances, which agrees with the ANN result and indicates local second-order clustering. At broader distances the line drops below the envelope, suggesting separation among some of the larger hotspots after the intensity surface has been controlled.

Interpretation and hypothesis

Overall, the wildfire pattern reflects a combination of first-order and second-order processes. The quadrat test and KDE show strong global variation in intensity, so environmental and human-access conditions differ across the study area. The ANN and short-distance L-function results show additional clustering among nearby events even after that global trend is included in the null model.

A plausible hypothesis is that wildfire ignitions cluster near roads, communities, recreation sites, and other accessible forest edges where human activity is concentrated. The North Carolina Forest Service reports that nearly 99% of the state’s wildfires are caused by human activities such as backyard burning and campfires. At the same time, the Croatan landscape includes fire-adapted longleaf pine habitat; the U.S. Forest Service explains that prescribed burning is used there to improve longleaf pine habitat and reduce hazardous vegetation. Together, human ignition opportunities and spatial differences in fuels could create broad hotspots, while repeated ignitions at the same accessible locations could produce the tight local clusters.

One limitation is that this file combines fires from many years and contains many repeated coordinates. Some locations may represent recurring fires, while others may be generalized reporting locations. Therefore, the small-distance clustering should not automatically be interpreted as one fire directly causing another.