PPOL 6805 / DSAN 6750: GIS for Spatial Data Science
Fall 2026
Wednesday, October 28, 2026
✅ Last Week (Oct 21): Framework for Formalizing Point \leadsto Areal Hypotheses
| Repulsion | ← | Complete Spatial Randomness | → | Clustering |
|---|
 |
 |
 |
source("../dsan-globals/_globals.r")
source("../dsan-globals/_globals-sf.r")
set.seed(6805)
library(tidyverse) |> suppressPackageStartupMessages()
library(sf) |> suppressPackageStartupMessages()
library(spatstat) |> suppressPackageStartupMessages()
library(mdthemes) |> suppressPackageStartupMessages()
library(leaflet) |> suppressPackageStartupMessages()You’ve been hired as an archaeologist – congratulations! Your job: determine whether arrangement of trees formed:
grid_points <- sf::st_sample(grid_buffer_sf, size=rep(1,100))
grid_buffer_sf |> ggplot(aes(shape='Cell')) +
geom_sf() +
geom_sf(data=grid_points) +
scale_shape_manual("Shape", values=c('Cell'=19)) +
theme_classic(base_size=ha_base) +
theme(
legend.title = element_blank(),
# legend.text = element_text(size=18)
) +
transparent_bg()
age <- runif(npoints(pois_ppp))
pair_dists <- pairdist(pois_ppp)
close <- (pair_dists < r)
later <- outer(age, age, ">")
killed <- apply(close & later, 1, any)
killed_ppp <- pois_ppp[killed]
alive_ppp <- pois_ppp[!killed]
pois_window_sf <- pois_ppp |> sf::st_as_sf() |> filter(label=="window")
pois_killed_sf <- killed_ppp |> sf::st_as_sf() |> filter(label=="point")
pois_alive_sf <- alive_ppp |> sf::st_as_sf() |> filter(label=="point")
alive_buff_sf <- pois_alive_sf |> sf::st_buffer(r) |> sf::st_union() |> sf::st_intersection(pois_window_sf)
ggplot() +
geom_sf(data=pois_window_sf) +
geom_sf(data=alive_buff_sf, aes(color='Inhibition', shape='Inhibition'), linetype='dashed') +
geom_sf(data=pois_killed_sf, aes(color='Dead', shape='Dead'), size=2, stroke=2) +
geom_sf(data=pois_alive_sf, aes(color='Alive', shape='Alive'), size=1, stroke=1) +
scale_shape_manual(name=NULL, values=c("Alive"=19, "Dead"=4, 'Inhibition'=21), labels=c("Alive", "Dead", "Inhibition")) +
scale_color_manual(name=NULL, values=c("Alive"="black", "Dead"=cb_palette[1], "Inhibition"="black"), labels=c("Alive", "Dead", "Inhibition")) +
guides(shape=guide_legend(override.aes=list(fill = "white"))) +
theme_classic(base_size = hn_base) +
theme(plot.margin = unit(c(0,0,0,0), "cm")) +
transparent_bg()
library(tidyverse)
library(spatstat)
set.seed(6809)
N <- 60
r_core <- 0.05
obs_window <- square(1)
### Clustered data
clust_ppp <- rMatClust(
kappa=6,
scale=r_core,
mu=10
)
clust_sf <- clust_ppp |> sf::st_as_sf()
clust_plot <- clust_sf |>
ggplot() +
geom_sf(size=2) +
theme_classic(base_size=18) +
transparent_bg()
ggsave("images/clust_ppp.png", clust_plot, width=3, height=3)
# Intensity fn
clust_intensity <- density(clust_ppp, sigma = 0.1)
png("images/clust_intensity.png")
par(mar=c(1.5,0,0,2), las=2, oma=c(0,0,0,0), cex=2, bg=NA)
plot(clust_intensity, main=NULL)
contour(clust_intensity, add = TRUE)
dev.off()
### PCF
clust_pcf <- spatstat.explore::pcf(
clust_ppp, divisor="d",
r=seq(from=0.00, to=0.50, by=0.01)
)
clust_pcf_plot <- clust_pcf |> ggplot(aes(x=r, y=iso)) +
geom_hline(yintercept=1, linetype='dashed', linewidth=1) +
geom_area(color='black', fill=cb_palette[1], alpha=0.75) +
scale_x_continuous(breaks=seq(from=0.0, to=1.0, by=0.1)) +
labs(x="Distance", y="Density") +
theme_classic(base_size=14) +
transparent_bg()
ggsave("images/clust_pcf.png", clust_pcf_plot, width=3, height=3)| Original Data | First-Order | Second-Order |
|---|---|---|
| N = 60 Events | Events modeled individually \Rightarrow Intensity Function \lambda(\mathbf{s}) |
Events modeled pairwise \Rightarrow Pairwise-Corr Function \textrm{pcf}(\vec{h}) |
 |
 |
 |
sq_base <- 16
sq_psize <- 2.5
obs_window <- square(1)
r0 <- 0.2
sq_df <- tibble::tribble(
~x, ~y,
0.5-r0,0.5-r0,
0.5+r0,0.5+r0,
0.5-r0,0.5+r0,
0.5+r0,0.5-r0
)
sq_sf <- sf::st_as_sf(
sq_df,
coords = c("x","y")
)
sq_ppp <- as.ppp(sq_sf, W=obs_window)
sq_ppp |> sf::st_as_sf() |> ggplot() +
geom_sf(size=sq_psize) +
theme_classic(base_size=sq_base) +
transparent_bg()
### PCF
pcf_result <- spatstat.explore::pcf(
sq_ppp,
divisor="d",
r=seq(from=0.00, to=0.8, by=0.01)
)
pcf_result |> ggplot(aes(x=r, y=iso)) +
geom_hline(yintercept=1, linetype='dashed', linewidth=1.5) +
geom_area(color='black', fill=cb_palette[1], alpha=0.75) +
scale_x_continuous(breaks=seq(from=0.0, to=1.0, by=0.1)) +
theme_classic(base_size=sq_base) +
transparent_bg()
csr_pcf_result <- spatstat.explore::pcf(
csr_ppp,
divisor="d",
r=seq(from=0.00, to=0.8, by=0.01)
)
csr_pcf_result |> ggplot(aes(x=r, y=iso)) +
geom_hline(yintercept=1, linetype='dashed', linewidth=1.5) +
geom_area(color='black', fill=cb_palette[1], alpha=0.75) +
scale_x_continuous(breaks=seq(from=0.0, to=1.0, by=0.1)) +
theme_classic(base_size=sq_base) +
transparent_bg()
sim_base <- 20
sim_psize <- 2
sim_xticks <- seq(from=0.0, to=1.0, by=0.2)
sim_yticks <- seq(from=0.0, to=1.0, by=0.2)
gen_pois_df <- function(num_sims=1) {
pois_sims <- spatstat.random::rpoispp(
lambda = 60, nsim=num_sims
)
return(tibble::as_tibble(pois_sims))
}
#pois_dfs <- gen_pois_df()
#pois_dfs |> head()
pois_sims <- spatstat.random::rpoispp(
lambda = 60, nsim=3
)
to_sim_df <- function(cur_sim, sim_name) {
cur_df <- tibble::as_tibble(cur_sim) |> mutate(sim=sim_name)
return(cur_df)
}
combined_df <- imap(.x=pois_sims, .f=to_sim_df) |> bind_rows()
combined_df |> ggplot(aes(x=x, y=y)) +
geom_point(size=sim_psize) +
facet_wrap(vars(sim)) +
coord_equal() +
theme_classic(base_size=sim_base) +
theme(panel.spacing.x = unit(2, "lines")) +
scale_x_continuous(breaks=sim_xticks) +
scale_y_continuous(breaks=sim_yticks) +
transparent_bg()
done:
r units of any existing point in \mathbf{S}
done=TRUE if \mathbf{S} has n points OR has been the same for giveup stepscapture.output(ssi_sims <- spatstat.random::rSSI(
r = 0.05, n=60, nsim=3
), file=nullfile())
combined_df <- imap(.x=ssi_sims, .f=to_sim_df) |> bind_rows()
combined_df |> ggplot(aes(x=x, y=y)) +
geom_point(size=sim_psize) +
facet_wrap(vars(sim)) +
coord_equal() +
theme_classic(base_size=24) +
theme(panel.spacing.x = unit(3, "lines")) +
scale_x_continuous(breaks=sim_xticks) +
scale_y_continuous(breaks=sim_yticks) +
transparent_bg()
scale centered at \mathbf{s}_imatclust_sims <- rMatClust(
kappa = 6,
scale = 0.075,
mu = 10,
nsim = 3
)
matclust_df <- imap(.x=matclust_sims, .f=to_sim_df) |> bind_rows()
matclust_plot <- matclust_df |> ggplot(aes(x=x, y=y)) +
geom_point(size=sim_psize) +
facet_wrap(vars(sim), nrow=1) +
coord_equal() +
theme_classic(base_size=sim_base) +
theme(panel.spacing.x = unit(2, "lines")) +
scale_x_continuous(breaks=sim_xticks) +
scale_y_continuous(breaks=sim_yticks) +
transparent_bg()
matclust_plotrMaternI() [Docs] |
rMaternII() [Docs] |
|---|---|
|
|
|
|
mI_sims <- rMaternI(
kappa = 60, r = 0.075, nsim=2
)
mII_sims <- rMaternII(
kappa = 60, r = 0.075, nsim=2
)
mI_combined_df <- imap(.x=mI_sims, .f=to_sim_df) |> bind_rows() |> mutate(sim=paste0("MI ",sim))
mII_combined_df <- imap(.x=mII_sims, .f=to_sim_df) |> bind_rows() |> mutate(sim=paste0("MII ",sim))
m_combined_df <- bind_rows(mI_combined_df, mII_combined_df)
m_plot <- m_combined_df |> ggplot(aes(x=x, y=y)) +
geom_point(size=sim_psize) +
facet_wrap(vars(sim), nrow=1) +
coord_equal() +
theme_classic(base_size=sim_base) +
theme(panel.spacing.x = unit(2, "lines")) +
scale_x_continuous(breaks=sim_xticks) +
scale_y_continuous(breaks=sim_yticks) +
transparent_bg()
m_plot
cox_pcol <- "black"
cox_bg <- "grey90"
cox_pch <- 21
# inhomogeneous LGCP with Gaussian covariance function
m <- as.im(function(x, y){
5 - 1.5 * (x - 0.5)^2 + 2 * (y - 0.5)^2
}, W=owin())
lgcp_sims <- rLGCP("gauss", m, var=0.15, scale =0.5, nsim=3)
# lgcp_combined_df <- imap(.x=ssi_sims, .f=to_sim_df) |> bind_rows()
plot_lgcp <- function(lgcp_sim) {
plot(attr(lgcp_sim, "Lambda"), main=NULL)
points(lgcp_sim, col=cox_pcol, bg=cox_bg, pch=cox_pch)
}
par(mfrow=c(1,3), mar=c(0,0,0,2), oma=c(0,0,0,0), las=2)
nulls <- lapply(X=lgcp_sims, FUN=plot_lgcp)
density() from spatstat makes “smart” choices for this bandwidth, so for now we just average the bandwidths estimated for cases and controls:[1] 11.45837To account for N_0 \neq N, normalize numerator and denominator:
\hat{\rho}(\mathbf{s}) = \frac{\hat{\lambda}(\mathbf{s}) / N}{\hat{\lambda}_0(\mathbf{s}) / N_0}
\Rightarrow \hat{\rho}(\mathbf{s}) = \frac{\hat{\lambda}(\mathbf{s})}{\alpha\hat{\lambda}_0(\mathbf{s})}, where \alpha = \frac{N}{N_0}
library(fields)
library(scico)
x <- int_cases$xcol
y <- int_cases$yrow
rr <- t(int_cases$v)/t(alpha * int_controls$v)
center <- 1.0
lim <- max(abs(rr - center), na.rm = TRUE)
# pal <- colorRampPalette(c("#2166AC", "white", "#B2182B"))(101)
par(
mar=c(6.5,3,1.5,0.5),
oma=c(1,1,1,1),
xpd=NA,
mgp = c(1.6, 0.6, 0)
)
pal <- scico(101, palette = "vik")
image.plot(x, y, rr, asp = 1, col=pal, zlim = c(center - lim, center + lim), main="Relative Risk Surface", xlab="Easting", ylab="Northing", bty='l')
HHHHHTHHHH…H, or could be fair (lots of heads by chance!)
Figure from Gimond (2024)
plot_rlabel <- function(rl_ppp) {
rl_win_sf <- window_to_sf(rl_ppp)
rl_points_sf <- points_to_sf(rl_ppp) |>
rename(
type = spatstat.geom..marks.x.
)
rl_plot <- ggplot() +
geom_sf(
data=rl_win_sf,
alpha=0.3
) +
geom_sf(
mapping=aes(fill=type, alpha=type),
data=rl_points_sf,
size=2,
# alpha=0.75,
shape=21,
stroke=0.1,
# color='white',
) +
# scale_color_manual(
# "Point Type",
# values=c("control"="gray","case"=cb_palette[1])
# ) +
scale_fill_manual(
"Type",
values=c("control"="gray","case"=cb_palette[1])
) +
scale_alpha_manual(
"Type",
values=c("control"=0.5, "case"=1.0)
) +
theme_classic(base_size=16) +
labs(
title="NE Britain Disease Cases"
) +
transparent_bg()
return(rl_plot)
}
set.seed(6810)
rl1_ppp <- rlabel(pbc)
plot_rlabel(rl1_ppp)
pbc_win_sf <- pbc |> sf::st_as_sf() |> filter(label == "window")
plot_pbc_sim <- function() {
pbc_sim <- spatstat.random::rpoispp(
lambda = alpha * int_controls
)
# Separate window and points
pbc_sf <- pbc_sim |> sf::st_as_sf()
pbc_point_sf <- pbc_sf |> filter(label == "point")
pbc_win_sf <- pbc_sf |> filter(label == "window")
pbc_plot <- ggplot() +
geom_sf(data=pbc_win_sf) +
geom_sf(data=pbc_point_sf, size=0.5, alpha=0.5) +
theme_classic() +
transparent_bg()
return(pbc_plot)
}
plot_pbc_sim()
ppp \leftrightarrow sf Conversionppp to sf Conversion:
| label | geom |
|---|---|
| window | POLYGON ((0 0, 1 0, 1 1, 0 … |
| point | POINT (0.7185641 0.5606081) |
| point | POINT (0.978131 0.9523815) |
| point | POINT (0.1744719 0.06257196) |
library(mdthemes) |> suppressPackageStartupMessages()
pois_sf |> ggplot() +
geom_sf(data=pois_sf |> filter(label=="window"), aes(fill='grey')) +
geom_sf(data=pois_sf |> filter(label != "window"), aes(color='black')) +
md_theme_classic(base_size=26) +
scale_fill_manual(name=NULL, values=c("gray90"), labels=c("<span style='font-family: mono'>label == 'window'</span>")) +
scale_color_manual(name=NULL, values=c("black"), labels=c("<span style='font-family: mono'>label == 'point'</span>")) +
labs(title = "<span style='font-family: mono'>ppp</span> → <span style='font-family: mono'>sf</span> Result") +
transparent_bg()
sf to ppp Conversion:
Planar point pattern: 100 points
window: polygonal boundary
enclosing rectangle: [0, 1] x [0, 1] unitstree_ppp_sf <- tree_ppp |> sf::st_as_sf()
tree_ppp_sf |> ggplot() +
geom_sf(aes(fill='gray90')) +
geom_sf(data=tree_ppp_sf |> filter(label != "window"), aes(color='black')) +
md_theme_classic(base_size=26) +
scale_fill_manual(name=NULL, values=c("gray90"), labels=c("<span style='font-family: mono'>tree_ppp$window</span>")) +
scale_color_manual(name=NULL, values=c("black"), labels=c("<span style='font-family: mono'>tree_ppp${x,y}</span>")) +
labs(
title = "<span style='font-family: mono'>sf</span> → <span style='font-family: mono'>ppp</span> Result",
x="<span style='font-family: mono'>tree_ppp$x</span>",
y="<span style='font-family: mono'>tree_ppp$y</span>"
) +
guides(
fill = guide_legend(order = 1),
color = guide_legend(order = 2)
) +
transparent_bg()