Code
source("../dsan-globals/_globals.r")PPOL 6805 / DSAN 6750: GIS for Spatial Data Science
Fall 2025
sf → ppp → spatial hypothesis testing!source("../dsan-globals/_globals.r")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)
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(0,0,0,2), las=2, oma=c(0,0,0,0), cex=2)
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)
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)
par(mar=c(0,0,0,2), las=2, oma=c(0,0,1,0))
sq_density <- density(sq_ppp)
plot(sq_density, main=NULL, xaxs='i', yaxs='i')
contour(sq_density, xaxs='i', yaxs='i', add = TRUE)
### 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)
csr_ppp <- spatstat.random::rpoispp(
lambda = 60,
win=obs_window
)
csr_ppp |> sf::st_as_sf() |> ggplot() +
geom_sf(size=sq_psize) +
theme_classic(base_size=sq_base)
par(
mar=c(0,0,0,2),
las=2,
oma=c(0,0,1,0)
)
csr_density <- density(csr_ppp)
plot(csr_density, main=NULL, xaxs="i", yaxs="i")
contour(csr_density, xaxs="i", yaxs="i", add = TRUE, lwd=1.5)
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)
ppp Objectsppp and sfppp Objectsspatstat book p. 41]sf objects, which contain data+geometries for any desired collection of \(N\) entities, ppp objects are required to have at least an observation window!sf Creation:
tree_df <- tibble::tibble(lon=runif(100,0,1), lat=runif(100,0,1), age=runif(100,0,1))
tree_sf <- sf::st_as_sf(
tree_df,
coords = c('lon', 'lat')
)
tree_sf |> head(4)| age | geometry |
|---|---|
| 0.8593935 | POINT (0.5035446 0.05519595) |
| 0.9770798 | POINT (0.2606818 0.3874983) |
| 0.9517643 | POINT (0.5464588 0.8271631) |
| 0.0592741 | POINT (0.7999682 0.5852131) |
ppp Creation:
pois_ppp <- spatstat.random::rpoispp(
lambda=100, win=spatstat.geom::square(1)
)
pois_pppPlanar point pattern: 93 points
window: rectangle = [0, 1] x [0, 1] units
attributes(pois_ppp)$names[1] "window" "n" "x" "y" "markformat"
pois_ppp$x |> head(4)[1] 0.51779807 0.62213583 0.84467836 0.05354204
ppp \(\leftrightarrow\) sf Conversionppp to sf Conversion:
pois_sf <- pois_ppp |> sf::st_as_sf()
pois_sf |> head(4)| label | geom |
|---|---|
| window | POLYGON ((0 0, 1 0, 1 1, 0 … |
| point | POINT (0.5177981 0.553718) |
| point | POINT (0.6221358 0.3933781) |
| point | POINT (0.8446784 0.8187611) |
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")
sf to ppp Conversion:
square_sfc <- sf::st_polygon(list(
matrix(c(0,0,1,0,1,1,0,1,0,0), nrow=5, byrow=TRUE)
)) |> sf::st_sfc()
tree_ppp <- as.ppp(
sf::st_as_sfc(tree_sf),
W=as.owin(square_sfc)
)
tree_pppPlanar point pattern: 100 points
window: polygonal boundary
enclosing rectangle: [0, 1] x [0, 1] units
tree_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))
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)
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(0,0,0,2), las=2, oma=c(0,0,0,0), cex=2)
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)
ggsave("images/clust_pcf.png", clust_pcf_plot, width=3, height=3)| Original Data | First-Order | Second-Order |
|---|---|---|
| \(N = 60\) Events | Events modeled individually \(\implies\) Intensity function \(\lambda(\mathbf{s})\) |
Events modeled pairwise \(\implies\) \(K\)-function \(K(\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)
par(mar=c(0,0,0,2), las=2, oma=c(0,0,1,0))
sq_density <- density(sq_ppp)
plot(sq_density, main=NULL, xaxs='i', yaxs='i')
contour(sq_density, xaxs='i', yaxs='i', add = TRUE)
### 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)
csr_ppp <- spatstat.random::rpoispp(
lambda = 60,
win=obs_window
)
csr_ppp |> sf::st_as_sf() |> ggplot() +
geom_sf(size=sq_psize) +
theme_classic(base_size=sq_base)
par(
mar=c(0,0,0,2),
las=2,
oma=c(0,0,1,0)
)
csr_density <- density(csr_ppp)
plot(csr_density, main=NULL, xaxs="i", yaxs="i")
contour(csr_density, xaxs="i", yaxs="i", add = TRUE, lwd=1.5)
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)