Week 6: Random Fields and Spatial Autocorrelation

PPOL 6805 / DSAN 6750: GIS for Spatial Data Science
Fall 2024

Author

Affiliation

Jeff Jacobs

jj1088@georgetown.edu

Open slides in new tab →

Random Fields

• (Now on slides rather than scribbled on the whiteboard!)

$D$ = Rectangle; $Z(\mathbf{s})$ here: height at coordinate $\mathbf{s}\in D$

Key Notation / Definition

• $d$-Dimensional “Spatial Process” (Schabenberger and Gotway 2004, 6)

$$\text{Data}=\{Z(\mathbf{s})\mid \mathbf{s}\in D\subset {\mathbb{R}}^d\}$$

• $d>1$: Data forms a Random Field (this class: $d=2$!)

	Geostatistical Data	Lattice/Region Data	Point Pattern
Criteria	Fixed $D$, Continuous	Fixed $D$, Discrete	Random subset $D^{\ast }\subseteq D$
Interest	Infer non-observed parts of $D$	Autocorrelation, clustering	Point-generating process
Example	$N$ trees ${\mathbf{s}}_1,{\mathbf{s}}_2,\dots ,{\mathbf{s}}_N$, observed within a sample window $D\subset {\mathbb{R}}^2$, ($D$ some finite plot of land) $Z({\mathbf{s}}_i)$: Attribute(s) at site ${\mathbf{s}}_i$ Example: Height. $Z({\mathbf{s}}_1)=500\text{m}$, $Z({\mathbf{s}}_2)=850\text{m}$, $\dots$	$Z(\mathbf{s})$ observed over $N\times N$ grid of plots $\Rightarrow$ Contiguity, Neighbors (next section of slides!) $\Rightarrow$ Autocorrelation: Are points around ${\mathbf{s}}_i$ likely to have values similar to $Z({\mathbf{s}}_i)$?	Unknown number of lightning strikes ${\mathbf{s}}_1,{\mathbf{s}}_2,\dots$ Contrast with geostatistical: all of $D$ is observed, but what determines the subset $D^{\ast }$ where events occur? “Unmarked”: Just locations “Marked”: Locations+info (e.g., intensity of strike)

Geostatistical Data: Dumping Sand onto a Table

• Domain $D$: tabletop; Random process: dumping sand onto the tabletop
• $\mathrm{\Omega }$: All possible realizations of process. ${\omega }_i\in \mathrm{\Omega }$ Particular realization

Figure 1: ${\omega }_1\in \mathrm{\Omega }$

Figure 2: ${\omega }_2\in \mathrm{\Omega }$

Figure 3: ${\omega }_3\in \mathrm{\Omega }$

Figure 4: ${\omega }_4\in \mathrm{\Omega }$

Geostatistical Modeling: Terrains Are Not Spiky Chaotic Landscapes!

Geostatistical Data $\to$ Point Pattern

• Indicator function: $\mathbb{1}[Z(\mathbf{s})>860\text{m}]$

set.seed(6805)
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(sf)

Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

library(terra)

terra 1.7.78

Attaching package: 'terra'

The following object is masked from 'package:tidyr':

    extract

library(latex2exp)
#### Define area
lat_range <- c(-5.0, -6.0)
lon_range <- c(15, 20)
#lon_center <- -5.4
#lat_center <- 36.75
# And the window around this centroid
lat_radius <- 0.05
lon_radius <- 0.1

gen_random_table <- function(world_num=1, return_data=FALSE, rand_seed=NULL) {
  if (!is.null(rand_seed)) {
    set.seed(rand_seed)
  } else {
    rand_seed <- sample(1:9999, size=1)
    set.seed(rand_seed)
  }
  lat_center <- runif(1, min=min(lat_range), max=max(lat_range))
  lon_center <- runif(1, min=min(lon_range), max=max(lon_range))
  lon_lower <- lon_center - lon_radius
  lon_upper <- lon_center + lon_radius
  lat_lower <- lat_center - lat_radius
  lat_upper <- lat_center + lat_radius
  coords <- data.frame(x = c(lon_lower, lon_lower, lon_upper, lon_upper),
                       y = c(lat_lower, lat_upper, lat_lower, lat_upper))
  coords_sf <- st_as_sf(coords, coords = c("x", "y"), crs = 4326)
  # Using geodata
  elev <- geodata::elevation_3s(lon = lon_center, lat = lat_center, res=2.5, path = getwd())
  elev <- terra::crop(elev, coords_sf)
  return(elev)
}

compute_hillshade <- function(elev) {
  # Calculate hillshade
  slopes <- terra::terrain(elev, "slope", unit = "radians")
  aspect <- terra::terrain(elev, "aspect", unit = "radians")
  hs <- terra::shade(slopes, aspect)
  return(hs)
}

plot_hillshade <- function(elev, hs, title=NULL) {
  ## Plot hillshading as basemap
  # (here using terra::plot, but could use tmap)
  # lat_upper_trimmed <- lat_upper - 0.001*lat_upper
  # lat_lower_trimmed <- lat_lower - 0.001*lat_lower
  base_plot <- terra::plot(
    hs, col = gray(0:100 / 100), legend = FALSE, axes = FALSE, mar=c(0,0,1,0), grid=TRUE
  )
  #    xlim = c(elev_xmin, elev_xmax), ylim = c(elev_ymin, elev_ymax))
  # overlay with elevation
  
  color_vec <- terrain.colors(25)
  plot(elev, col = color_vec, alpha = 0.5, legend = FALSE, axes = FALSE, add = TRUE)
  
  # add contour lines
  terra::contour(elev, col = "grey30", add = TRUE)
  if (!is.null(title)) {
    world_name <- TeX(sprintf(r"(World $\omega_{%d} \in \Omega$)", world_num))
    title(world_name)
  }
}
elev_r1 <- gen_random_table(rand_seed=6810)
hs_r1 <- compute_hillshade(elev_r1)
plot_hillshade(elev_r1, hs_r1)

elev_indic_r1 <- terra::ifel(
  elev_r1 > 860, elev_r1, NA
)
plot(elev_indic_r1)

elev_r2 <- gen_random_table(rand_seed=6815)
hs_r2 <- compute_hillshade(elev_r2)
plot_hillshade(elev_r2, hs_r2)

my_quant <- function(x) quantile(x, 0.99)
elev_99 <- terra::global(elev_r2, my_quant)[1,1]
elev_indic_r2 <- terra::ifel(
  elev_r2 > elev_99, elev_r2, NA
)
plot(elev_indic_r2)

Geospatial Data $\to$ Lattice/Region Data

table_elev <- gen_random_table(rand_seed=6805)
table_hs <- compute_hillshade(table_elev)
#plot_hillshade(table_elev, table_hs)
#par(mfrow=c(1,2), mar=c(0,0,1,1))
plot(table_elev)

(Pretend this is infinite-resolution) Infinite number of neighbors around each point!

elev_coarse <- terra::aggregate(table_elev, 12)
plot(elev_coarse)

Finite number of neighbors around each observation $⟹$ Can study correlations between neighboring cells!

Doing Things with nbdep

• “Neighbor” Definitions
• Testing Hypotheses About Neighborhoods!

Who Are My Neighbors?

• Contiguity-Based:

• Distance-Based

• $K$-Nearest Neighbors

Distance-Based Neighbors

From Xie (2023)

Spatial Autocorrelation

• Location $i$ has high value $⟹$ locations near $i$ more likely to have high values

Moran’s $I$

$$I=\underset{\text{Inverse of Variance}}{\boxed{{\displaystyle \frac{n}{\sum _{i=1}^n(y_i-\bar{y}{)}^2}}}}\frac{\stackrel{\text{Weighted Covariance}}{\overbrace{\sum _{i=1}^n\sum _{j=1}^n{w}_{ij}(y_i-\bar{y})(y_j-\bar{y})}}}{\underset{\text{Normalize Weights}}{\underbrace{\sum _{i=1}^n\sum _{j=1}^n{w}_{ij}}}}$$

• $I$ is Large when:
  • $y_i$ and $y_j$ are neighbors: $w_{ij}$, and
  • $y_i$ and $y_j$ large at the same time: $(y_i-\bar{y})(y_j-\bar{y})$

Local Indicators of Spatial Autocorrelation (LISA)

• tldr: See how much a given location $\mathbf{s}$ contributes to overall Moran’s $I$
• Local Moran’s $I$:

$$I_i=\frac{y_i-\bar{y}}{S_i^2}\sum _{j=1}^n{w}_{ij}(y_j-\bar{y})$$

Our Null Hypothesis: Complete Spatial Randomness (CSR)

Leaving You With a Challenge

• These point patterns can be classified as clustered or regular:

Leaving You With a Challenge

References

Schabenberger, Oliver, and Carol A. Gotway. 2004. Statistical Methods for Spatial Data Analysis. CRC Press.

Xie, Sherrie. 2023. “Analyzing Geospatial Data in R.”