Tutorial

Submit your knitted .html to Canvas. Completion of the mini-challenge(s) is graded on the meaningful attempt rubric (1 point)

Some phenomena are spatially continuous, meaning they are present at every location rather than confined to discrete points.

In many cases, however, data about these continuous variables are only collected at discrete locations (i.e. weather monitoring stations, stream gauges, sample sites). Therefore, an important task is to predict these continuous values at unsampled locations. This is called spatial interpolation. In this tutorial we will explore several strategies for interpolating

In this tutorial, we will explore several strategies for interpolating data. Ultimately, the choice of interpolation method depends on several factors, including the strength and scale of autocorrelation, the density of sampling locations, the need for estimates of uncertainty, and the computational complexity.

We will use data from weather stations across North Carolina during a 2025 heat wave (July 30, 2025). The data is from the North Carolina State Climate Office

Start by creating a new copy of the class .Rmd template and saving it as “continuous_surfaces.Rmd” As you move through the tutorial, copy and paste the commands into your own .Rmd file. Any time you see a question in the document these are meant for personal reflection (you won’t turn these answers in)

Reading in Data

library(sf)
library(tidyverse)
library(gstat)
library(sp)
library(terra)
library(tmap)
library(INLA)

options(scipen = 999)

#heatwave data
heatwave_temp <- st_read("https://drive.google.com/uc?export=download&id=1tZQU5NwZ3Z1lBtFbJTvM7clB1z0sXMI9") |> distinct()

#nc boundaries
nc_boundaries <- st_read("https://drive.google.com/uc?export=download&id=1SNy838Q7f7rvL6TFLmIZm-sKyRO7Htjw")

Spatial Structure

To predict continuous surfaces, we must assume that the underlying phenomenon has a predictable spatial structure. Methods for surface prediction rely on spatial autocorrelation. When a variable shows little or no spatial autocorrelation, it is not a good candidate for continuous surface prediction.

Let’s first look at our temperature values:

tm_shape(heatwave_temp) + tm_dots("value", tm_scale_intervals(style = "quantile", values = "-brewer.rd_yl_gn")) + tm_shape(nc_boundaries) + tm_borders()

Just based on the pattern visible on the map, this variable clearly is autocorrelated. However, if we were unsure, we could calculate a Moran’s I.

Methods for Interpolation

IDW

Inverse Distance Weighting (IDW) is a deterministic interpolation method that predicts values at unsampled locations as a weighted average of nearby observations. Points closer to the prediction location receive more weight than points farther away. IDW is a simple method that is often used to generate an exploratory surface.

IDW does not estimate an underlying statistical model of spatial dependence. Instead, it assumes that similarity decreases with distance and uses that assumption directly to compute predictions.

The main parameter in IDW is the inverse distance power (IDP), which controls how rapidly influence declines as distance increases. Lower IDP (closer to 1) values allow more distant observations to contribute to predictions, while higher IDP (greater than 3) values emphasize local variation and weight nearby points more heavily.

To select an appropriate IDP, we can compare several reasonable values using something called k-fold cross-validation. This approach iteratively divides the dataset into training and testing subsets and withholds part of the data. Predictions are generated using the remaining observations and value with the lowest prediction error (usually using RMSE) is selected.

#turn sf objects into sp objects for gstats
heatwave_sp <- as(heatwave_temp, "Spatial")
nc_sp <- as(nc_boundaries, "Spatial")

#test values
idp_vals <- c(1, 1.5, 2, 2.5, 3, 4)

#store results
idw_cv_results <- data.frame(
  idp = idp_vals,
  ME = NA,
  MAE = NA,
  RMSE = NA
)

#loop that calculates prediction error at all of the test vallues
for (i in seq_along(idp_vals)) {
  
  idw_model <- gstat(
    formula = value ~ 1, #idw formula
    locations = heatwave_sp,
    set = list(idp = idp_vals[i])
  )
  
  cv <- gstat.cv(
    object = idw_model,
    nfold = 5   #5-fold cross-validation
  )
  
  idw_cv_results$ME[i] <- mean(cv$residual, na.rm = TRUE)
  idw_cv_results$MAE[i] <- mean(abs(cv$residual), na.rm = TRUE)
  idw_cv_results$RMSE[i] <- sqrt(mean(cv$residual^2, na.rm = TRUE))
}

#view benchmarking results
idw_cv_results

##   idp        ME      MAE     RMSE
## 1 1.0 0.2013001 2.219834 3.109222
## 2 1.5 0.2971721 1.973363 2.878005
## 3 2.0 0.3207913 1.944057 2.777403
## 4 2.5 0.2828201 2.089997 2.998719
## 5 3.0 0.1712999 2.106596 3.103021
## 6 4.0 0.2149711 2.231425 3.165803

Q1: Based on our benchmarking results, what IDP should we choose? (Note that the RMSE might vary slightly depending on the simulations)

Answer: I would choose an IDP of 2 because it has the lowest RMSE (2.777) and the lowest MAE (1.944) in these cross-validation results. Its mean error is fairly close to zero as well. Since the folds are random, the exact values can change slightly between runs, but IDP = 2 is the best choice for the results shown here.

To actually create a prediction surface, we first need to create a prediction grid, which defines the locations where we want to estimate temperature values. The grid (raster) divides the study area into cells, and the interpolation method then will generate a predicted value for each cell. The size of each pixel determines the spatial resolution of the final surface. Smaller pixels produce a smoother, more detailed surface but require more computation, while larger pixels produce a coarser surface with less spatial detail. In this tutorial, we create a grid that covers the North Carolina boundary and will use this same grid for all interpolation methods so that the resulting prediction surfaces can be compared directly.

# create a regular prediction grid. We say we want 2000 pixels
grd <- as.data.frame(spsample(nc_sp, "regular", n = 2000))

#set grid coordinates
names(grd) <- c("X", "Y")
coordinates(grd) <- c("X", "Y")
gridded(grd) <- TRUE
fullgrid(grd) <- TRUE

#assign CRS (which is 2264, you MUST use a projected coordinate system for interpolation)
crs(grd) <- crs(heatwave_sp)

# run IDW interpolation using the best IDP from cross-validation
heatwave_idw <- idw(
  formula = value ~ 1,
  locations = heatwave_sp,
  newdata = grd,
  idp = 2
)

## [inverse distance weighted interpolation]

#rasterize interpolation output
r <- rast(heatwave_idw)

#mask to North Carolina boundaries
r_m <- mask(r, nc_boundaries)

#map interpolated temperatures
tm_shape(r_m) +
  tm_raster(
    col = "var1.pred",
    col.scale = tm_scale_intervals(
      values = "-rd_yl_gn",
      style = "quantile"
    )
  )

Q2: Re-run the IDW using a different IDP and compare the results.

#compare the selected IDP of 2 with a more locally focused IDP of 4
heatwave_idw_4 <- idw(
  formula = value ~ 1,
  locations = heatwave_sp,
  newdata = grd,
  idp = 4
)

r_4 <- rast(heatwave_idw_4)
r_4_m <- mask(r_4, nc_boundaries)

tm_shape(r_4_m) +
  tm_raster(
    col = "var1.pred",
    col.scale = tm_scale_intervals(
      values = "-rd_yl_gn",
      style = "quantile"
    )
  )

Heat-wave IDW surface using IDP 4

Answer: With IDP = 4, nearby stations receive much more weight, so the map has sharper local hot and cool spots around individual stations. The IDP = 2 surface is smoother because observations farther away still influence each prediction. The higher-power version preserves more local variation, but it also creates more of a bull’s-eye pattern and had a higher cross-validation error than IDP = 2.

Kriging

Kriging assumes that the observed data are realizations of an underlying spatial random field, meaning the values we measure at specific locations are part of a continuous stochastic process across space. Unlike IDW, which determines prediction weights solely based on distance, kriging uses a variogram model to describe how strongly observations are related to one another across space. This model characterizes how spatial similarity declines as the distance between locations increases, and kriging uses this relationship to determine the weights applied to nearby observations when predicting values at unsampled locations.

To understand the spatial structure of the data, we first calculate an empirical variogram on our data. The empirical variogram summarizes how differences between observed values change with distance by computing the average semivariance between pairs of observations at different separation distances.

#compute empirical variogram
v <- variogram(value ~ 1, data = heatwave_sp)

#plot empirical variogram
plot(v)

We then match the overall shape of the empirical variogram to a theoretical variogram model, which provides a smooth mathematical function describing the spatial dependence structure of the data. The most commonly used models include the spherical, exponential, and Gaussian variogram models.

show.vgms(par.strip.text = list(cex = 0.75))

In our case, we will use a spherical variogram model.

Variogram models are typically described using three main parameters: the nugget, the sill, and the range. When fitting a variogram model, we provide initial estimates for these parameters so the fitting algorithm has a starting point. The goal is to roughly match the overall shape of the empirical variogram to one of the available theoretical models.

Once these starting values are specified, the fitting procedure refines the parameters to produce a variogram model that best represents the spatial structure observed in the data. This fitted model is then used by the kriging algorithm to determine the weights applied to nearby observations when generating predictions.

In our case, we will define the nugget at 10, the sill at 17, and the range at 150000

#specify an initial variogram model
vgm_initial <- vgm(
  psill = 17,
  model = "Sph",
  range = 150000,
  nugget = 10
)

#fit empirical to initial model
vgm_fit <- fit.variogram(v, vgm_initial)

#plot fitted variogram model 
plot(v, vgm_fit)

#fit kriging
krig_pred <- krige(
  formula = value ~ 1,
  locations = heatwave_sp,
  newdata = grd,
  model = vgm_fit
)

## [using ordinary kriging]

#rasterize predictions
krig_rast <- rast(krig_pred)

#mask to North Carolina boundary
krig_mask <- mask(krig_rast, nc_boundaries)

#map kriging predictions
tm_shape(krig_mask) +
  tm_raster(
    col = "var1.pred",
    col.scale = tm_scale_intervals(
      values = "-rd_yl_gn",
      style = "quantile"
    )
  )

# map kriging prediction variance
tm_shape(krig_mask) +
  tm_raster(
    col = "var1.var",
    col.scale = tm_scale_intervals(
      values = "viridis",
      style = "quantile"
    )
  )

Q3. Compare our kriging results to our IDW. What looks different?

Answer: Both maps show the same broad pattern, including cooler temperatures in western North Carolina and warmer temperatures across much of the central part of the state. The IDW map has a large, smooth hot region but also shows small bull’s-eye patterns around some stations. The kriging map has more irregular pockets and transitions because its weights come from the fitted variogram rather than distance alone. Kriging also slightly reduces the most extreme predictions.

Q4. Look at the kriging variance. This represents the uncertainty of the prediction at that location. Taking the square of the variance is the standard error (i.e. if the variance is 17, that means the standard error would be about 4 degrees). What could cause differences in variance across the study region?

Answer: Prediction variance is generally lower near stations and in places where stations are densely clustered because the model has more nearby information. It becomes higher where stations are sparse, where a location is far from its nearest station, and near the state boundary because there are fewer observations surrounding the prediction location. The fitted nugget, sill, and range also affect the amount and spatial pattern of uncertainty. Technically, the standard error is the square root of the variance, so a variance of 17 corresponds to a standard error of about 4.1 degrees.

Mini Challenge

#mean maximum jan temp
jan_temp <- st_read("https://drive.google.com/uc?export=download&id=17abGtjWgvVWw8Wv69lvT_Q9qbUA9eblm")

## Reading layer `mean_maximum_jan_temp' from data source 
##   `https://drive.google.com/uc?export=download&id=17abGtjWgvVWw8Wv69lvT_Q9qbUA9eblm' 
##   using driver `GeoJSON'
## Simple feature collection with 42 features and 1 field
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 946463.5 ymin: 191099.5 xmax: 3008377 ymax: 989924.7
## Projected CRS: NAD83 / North Carolina (ftUS)

Apply IDW and Kriging interpolation approaches to estimate mean maximum January daily temperature at unmeasured locations across North Carolina. Make sure to use an appropriate IDP for your IDW and select an appropriate variogram model and nugget, sill, and range estimations.

IDW interpolation

I first compared several IDP values using five-fold cross-validation. I used a fixed random seed so that the comparison is reproducible.

jan_temp_sp <- as(jan_temp, "Spatial")

set.seed(734)
jan_idp_vals <- c(1, 1.5, 2, 2.5, 3, 4)
jan_idw_cv <- data.frame(
  idp = jan_idp_vals,
  ME = NA_real_,
  MAE = NA_real_,
  RMSE = NA_real_
)

for (i in seq_along(jan_idp_vals)) {
  jan_idw_model <- gstat(
    formula = value ~ 1,
    locations = jan_temp_sp,
    set = list(idp = jan_idp_vals[i])
  )

  jan_cv <- gstat.cv(jan_idw_model, nfold = 5)
  jan_idw_cv$ME[i] <- mean(jan_cv$residual, na.rm = TRUE)
  jan_idw_cv$MAE[i] <- mean(abs(jan_cv$residual), na.rm = TRUE)
  jan_idw_cv$RMSE[i] <- sqrt(mean(jan_cv$residual^2, na.rm = TRUE))
}

jan_idw_cv
best_idp <- jan_idw_cv$idp[which.min(jan_idw_cv$RMSE)]
best_idp

##   idp     ME   MAE  RMSE
## 1 1.0 -0.154 1.370 1.940
## 2 1.5 -0.186 1.194 1.772
## 3 2.0 -0.196 1.139 1.691
## 4 2.5 -0.191 1.113 1.675
## 5 3.0 -0.179 1.110 1.694
## 6 4.0 -0.145 1.148 1.769
## [1] 2.5

IDP = 2.5 produced the lowest RMSE, so I used it for the final January IDW surface.

jan_idw <- idw(
  formula = value ~ 1,
  locations = jan_temp_sp,
  newdata = grd,
  idp = best_idp
)

jan_idw_rast <- rast(jan_idw)
jan_idw_mask <- mask(jan_idw_rast, nc_boundaries)

tm_shape(jan_idw_mask) +
  tm_raster(
    col = "var1.pred",
    col.scale = tm_scale_intervals(
      values = "-rd_yl_gn",
      style = "quantile"
    )
  )

IDW interpolation of mean maximum January temperature

Kriging interpolation

The empirical variogram rises with distance and then levels off, so I selected a spherical model. Approximate starting values were a nugget of 0.35, partial sill of 4.22, and range of 745,000 feet. The fitting step then refines those starting estimates.

jan_v <- variogram(value ~ 1, data = jan_temp_sp)

jan_vgm_initial <- vgm(
  psill = 4.22,
  model = "Sph",
  range = 745000,
  nugget = 0.35
)

jan_vgm_fit <- fit.variogram(jan_v, jan_vgm_initial)
plot(jan_v, jan_vgm_fit)

Empirical January temperature variogram with fitted spherical model

jan_krig <- krige(
  formula = value ~ 1,
  locations = jan_temp_sp,
  newdata = grd,
  model = jan_vgm_fit
)

jan_krig_rast <- rast(jan_krig)
jan_krig_mask <- mask(jan_krig_rast, nc_boundaries)

tm_shape(jan_krig_mask) +
  tm_raster(
    col = "var1.pred",
    col.scale = tm_scale_intervals(
      values = "-rd_yl_gn",
      style = "quantile"
    )
  )

Ordinary kriging interpolation of mean maximum January temperature

tm_shape(jan_krig_mask) +
  tm_raster(
    col = "var1.var",
    col.scale = tm_scale_intervals(
      values = "viridis",
      style = "quantile"
    )
  )

Kriging prediction variance for mean maximum January temperature

Mini-challenge interpretation: Both methods show a clear west-to-east temperature pattern. The coolest mean January maximum temperatures are in the northwestern mountains, while the warmest values are concentrated across the southern Piedmont and Coastal Plain. IDW follows individual stations more closely and predicts a wider range, while kriging produces a smoother regional surface based on the fitted spatial structure. Kriging uncertainty is lowest near groups of stations and highest in the far western edge and other areas with limited nearby observations.

Tutorial: Continuous Surfaces