Spatial Interpolation of Point Data in R

Introduction

Spatial interpolation is a technique used to estimate values at unmeasured locations based on values from sampled locations. This is particularly useful in environmental sciences, geology, meteorology, and many other fields where collecting data at every location is impractical or impossible.

In this article, we’ll explore several interpolation methods available in R and demonstrate their implementation with practical examples.

Prerequisites

We’ll be using several R packages for spatial interpolation:

# Install required packages
install.packages(c("sf", "sp", "gstat", "stars", "ggplot2", "dplyr", "viridis"))

# Load libraries
library(sf)
library(sp)
library(gstat)
library(stars)
library(ggplot2)
library(dplyr)
library(viridis)

Common Interpolation Methods

1. Inverse Distance Weighting (IDW)

IDW is one of the simplest interpolation methods. It assumes that the influence of a sampled point decreases with distance, with closer points having more influence than distant ones.

Advantages: - Simple and intuitive - No statistical assumptions required - Fast computation

Disadvantages: - Can produce bull’s-eye patterns - No assessment of prediction uncertainty - Sensitive to clustering of sample points

Implementation

# Create sample point data
set.seed(123)
n <- 50
points_df <- data.frame(
  x = runif(n, 0, 100),
  y = runif(n, 0, 100),
  value = rnorm(n, 50, 10)
)

# Convert to spatial object
points_sf <- st_as_sf(points_df, coords = c("x", "y"), crs = NA)

# Create a grid for interpolation
grid <- st_bbox(points_sf) %>%
  st_as_stars(dx = 1, dy = 1) %>%
  st_set_crs(NA)

# Perform IDW interpolation
idw_result <- idw(value ~ 1, 
                  locations = as(points_sf, "Spatial"),
                  newdata = as(grid, "Spatial"),
                  idp = 2)  # power parameter

# Convert back to stars object
idw_stars <- st_as_stars(idw_result)

# Plot results
ggplot() +
  geom_stars(data = idw_stars, aes(fill = var1.pred, x = x, y = y)) +
  geom_sf(data = points_sf, size = 1, color = "black") +
  scale_fill_viridis(name = "Predicted\nValue") +
  theme_minimal() +
  labs(title = "IDW Interpolation Results",
       x = "X Coordinate", y = "Y Coordinate")

2. Kriging

Kriging is a geostatistical method that provides the best linear unbiased prediction. It uses variogram modeling to capture spatial autocorrelation and provides uncertainty estimates.

Types of Kriging: - Ordinary Kriging: Assumes constant but unknown mean - Universal Kriging: Accounts for trends in the data - Simple Kriging: Assumes known constant mean

Advantages: - Provides uncertainty estimates - Accounts for spatial autocorrelation - Generally more accurate than IDW

Disadvantages: - More complex - Requires variogram modeling - Computationally intensive for large datasets

Implementation

# Convert to sp objects for gstat
points_sp <- as(points_sf, "Spatial")
grid_sp <- as(grid, "Spatial")

# Compute empirical variogram
variogram_emp <- variogram(value ~ 1, points_sp)

# Fit variogram model
variogram_model <- fit.variogram(variogram_emp, 
                                 model = vgm(psill = 100, 
                                            model = "Sph", 
                                            range = 30, 
                                            nugget = 10))

# Plot variogram
plot(variogram_emp, variogram_model, 
     main = "Empirical and Fitted Variogram")

# Perform ordinary kriging
krige_result <- krige(value ~ 1,
                      locations = points_sp,
                      newdata = grid_sp,
                      model = variogram_model)

# Convert to stars and plot
krige_stars <- st_as_stars(krige_result)

ggplot() +
  geom_stars(data = krige_stars, aes(fill = var1.pred, x = x, y = y)) +
  geom_sf(data = points_sf, size = 1, color = "black") +
  scale_fill_viridis(name = "Predicted\nValue") +
  theme_minimal() +
  labs(title = "Ordinary Kriging Results",
       x = "X Coordinate", y = "Y Coordinate")

# Plot prediction variance
ggplot() +
  geom_stars(data = krige_stars, aes(fill = var1.var, x = x, y = y)) +
  scale_fill_viridis(name = "Prediction\nVariance") +
  theme_minimal() +
  labs(title = "Kriging Prediction Variance",
       x = "X Coordinate", y = "Y Coordinate")

3. Thin Plate Splines

Thin plate splines are a smoothing technique that minimizes the bending energy of the interpolating surface.

Advantages: - Produces smooth surfaces - No arbitrary parameters - Works well for scattered data

Disadvantages: - Can be computationally expensive - May oversmooth in some cases

Implementation

library(fields)

# Perform thin plate spline interpolation
tps_model <- Tps(x = points_df[, c("x", "y")], 
                 Y = points_df$value)

# Create prediction grid
grid_df <- expand.grid(
  x = seq(0, 100, length.out = 100),
  y = seq(0, 100, length.out = 100)
)

# Predict
grid_df$predicted <- predict(tps_model, 
                             x = grid_df[, c("x", "y")])

# Plot
ggplot() +
  geom_raster(data = grid_df, aes(x = x, y = y, fill = predicted)) +
  geom_point(data = points_df, aes(x = x, y = y), size = 1) +
  scale_fill_viridis(name = "Predicted\nValue") +
  theme_minimal() +
  labs(title = "Thin Plate Spline Interpolation",
       x = "X Coordinate", y = "Y Coordinate")

Comparing Methods

When choosing an interpolation method, consider:

1. Data characteristics: Are your data clustered or evenly distributed?
2. Spatial dependence: Is there strong spatial autocorrelation?
3. Uncertainty: Do you need uncertainty estimates?
4. Computational resources: How large is your dataset?
5. Smoothness: Do you expect smooth or abrupt changes?

Cross-Validation

Cross-validation helps compare interpolation methods:

# Leave-one-out cross-validation for IDW
idw_cv <- krige.cv(value ~ 1, 
                   locations = points_sp, 
                   nfold = nrow(points_sp))

# Calculate RMSE
idw_rmse <- sqrt(mean(idw_cv$residual^2))

# Cross-validation for kriging
krige_cv <- krige.cv(value ~ 1,
                     locations = points_sp,
                     model = variogram_model,
                     nfold = nrow(points_sp))

krige_rmse <- sqrt(mean(krige_cv$residual^2))

cat("IDW RMSE:", round(idw_rmse, 2), "\n")
cat("Kriging RMSE:", round(krige_rmse, 2), "\n")

Best Practices

1. Explore your data: Plot the sample locations and values before interpolation
2. Check for trends: Detrend data if necessary before kriging
3. Validate variogram models: Ensure the fitted variogram makes sense
4. Use cross-validation: Compare methods objectively
5. Consider anisotropy: Account for directional differences in spatial correlation
6. Check edge effects: Interpolation near boundaries is less reliable
7. Document parameters: Keep track of power parameters, variogram models, etc.

Conclusion

Spatial interpolation is a powerful tool for estimating values across a continuous surface from point samples. Each method has its strengths and appropriate use cases:

• Use IDW for quick estimates when you don’t need uncertainty measures
• Use Kriging when you need optimal predictions with uncertainty estimates
• Use Thin Plate Splines when you expect smooth surfaces

Always validate your interpolation results and choose methods appropriate for your data and research questions.

Additional Resources

• The gstat package documentation

---

📧 Contact & Collaboration

Have questions about this analysis or interested in collaborating on geospatial projects? We’d love to hear from you!

Get in touch with our research team: - Email: mapcrafty@gmail.com - Subject line: “Inquiry about Remote Sensing Applications for urban growth Monitoring”

Whether you’re working on similar research, need technical consultation, or want to discuss potential collaborations in geospatial analysis, don’t hesitate to reach out. Our team is always excited to connect with fellow researchers and practitioners in the GIS and remote sensing community.

We typically respond within 24-48 hours and welcome discussions about methodology, data sources, and potential research partnerships.