SPATIAL INTERPOLATION - PARAGUAY SOIL
MOISTURE DATA
0. LOAD PACKAGES
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.2.1 ✔ readr 2.2.0
## ✔ forcats 1.0.1 ✔ stringr 1.6.0
## ✔ ggplot2 4.0.2 ✔ tibble 3.3.1
## ✔ lubridate 1.9.5 ✔ tidyr 1.3.2
## ✔ purrr 1.2.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(sp)
library(gstat)
1. READ DATA
coordinates_data <- read.csv(
"~/Downloads/Manzano2022COMPAG_coordinates.csv"
)
measurements_data <- read.csv(
"~/Downloads/Manzano2022COMPAG_data.csv"
)
# Check data structure
str(coordinates_data)
## 'data.frame': 12 obs. of 3 variables:
## $ Node_ID: int 1 2 3 4 5 6 7 8 9 10 ...
## $ x : num 0 2.21 22.99 25.7 31.3 ...
## $ y : num 0 -9.16 3.9 -5.29 -11.43 ...
str(measurements_data)
## 'data.frame': 15777 obs. of 4 variables:
## $ t : chr "09-Mar-2022 13:02:32" "09-Mar-2022 13:09:24" "09-Mar-2022 13:09:51" "09-Mar-2022 13:11:50" ...
## $ Node_ID: int 9 3 3 10 1 9 4 4 3 5 ...
## $ Temp : num 27.1 27.4 27.4 27.4 27.1 ...
## $ Moist : num 41.2 38.4 38.4 42 56.7 ...
head(coordinates_data)
## Node_ID x y
## 1 1 0.000 0.000
## 2 2 2.210 -9.159
## 3 3 22.992 3.905
## 4 4 25.696 -5.292
## 5 5 31.297 -11.428
## 6 6 47.449 -5.017
head(measurements_data)
## t Node_ID Temp Moist
## 1 09-Mar-2022 13:02:32 9 27.06 41.22
## 2 09-Mar-2022 13:09:24 3 27.41 38.37
## 3 09-Mar-2022 13:09:51 3 27.41 38.37
## 4 09-Mar-2022 13:11:50 10 27.41 42.04
## 5 09-Mar-2022 13:19:28 1 27.06 56.73
## 6 09-Mar-2022 13:22:27 9 27.41 41.22
# Check number of observations and sensor locations
cat("Number of measurement rows:",
nrow(measurements_data), "\n")
## Number of measurement rows: 15777
cat("Number of unique nodes:",
length(unique(measurements_data$Node_ID)), "\n")
## Number of unique nodes: 12
2. DATA PREPROCESSING
# Aggregate repeated measurements into one mean per sensor node
# This creates one spatial observation per unique location.
mean_data <- measurements_data %>%
group_by(Node_ID) %>%
summarise(
Moist = mean(Moist, na.rm = TRUE),
Temp = mean(Temp, na.rm = TRUE),
n_obs = n(),
.groups = "drop"
)
# Merge with spatial coordinates
data_spatial <- mean_data %>%
left_join(coordinates_data, by = "Node_ID") %>%
select(Node_ID, x, y, Moist, Temp, n_obs)
# Remove missing values
data_spatial <- data_spatial %>%
filter(
!is.na(x),
!is.na(y),
!is.na(Moist)
)
# Check duplicated coordinates
duplicate_coordinates <- data_spatial %>%
count(x, y) %>%
filter(n > 1)
print(duplicate_coordinates)
## # A tibble: 0 × 3
## # ℹ 3 variables: x <dbl>, y <dbl>, n <int>
# Check final dataset
print(data_spatial)
## # A tibble: 12 × 6
## Node_ID x y Moist Temp n_obs
## <int> <dbl> <dbl> <dbl> <dbl> <int>
## 1 1 0 0 64.1 24.7 1319
## 2 2 2.21 -9.16 46.1 25.6 1737
## 3 3 23.0 3.90 47.3 25.2 1310
## 4 4 25.7 -5.29 64.3 25.3 1338
## 5 5 31.3 -11.4 45.7 24.5 1325
## 6 6 47.4 -5.02 43.9 24.8 1275
## 7 7 58.8 -5.86 44.1 24.9 1331
## 8 8 63.8 0.334 46.1 25.4 1203
## 9 9 74.0 0.019 48.5 25.3 1305
## 10 10 77.1 -10.9 52.6 25.7 1165
## 11 11 85.8 1.21 47.5 25.2 1184
## 12 12 90.3 -7.53 42.2 25.4 1285
cat("Final number of spatial observations:",
nrow(data_spatial), "\n")
## Final number of spatial observations: 12
3. EXPLORATORY DATA
ANALYSIS
# Summary statistics
summary(data_spatial$Moist)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 42.19 45.28 46.71 49.38 49.56 64.33
# Spatial distribution of soil moisture
ggplot(
data_spatial,
aes(x = x, y = y, color = Moist)
) +
geom_point(size = 4) +
geom_text(
aes(label = Node_ID),
nudge_y = 0.5,
size = 3
) +
scale_color_viridis_c() +
coord_equal() +
labs(
title = "Spatial Distribution of Soil Moisture",
subtitle = "Paraguay Sensor Network",
x = "X Coordinate (m)",
y = "Y Coordinate (m)",
color = "Mean Moisture"
) +
theme_minimal()
4. DEFINE INTERPOLATION
FUNCTIONS
# Function 1: Closest Observation
# Prediction = value of the nearest observed point
predict_closest <- function(train, test_point) {
distances <- sqrt(
(train$x - test_point$x)^2 +
(train$y - test_point$y)^2
)
nearest_index <- which.min(distances)
return(train$Moist[nearest_index])
}
# Function 2: Nearest Neighbors
# Prediction = average of k nearest observations
# Here k = 3
predict_nn <- function(train, test_point, k = 3) {
distances <- sqrt(
(train$x - test_point$x)^2 +
(train$y - test_point$y)^2
)
k <- min(k, nrow(train))
nearest_indices <- order(distances)[1:k]
prediction <- mean(
train$Moist[nearest_indices]
)
return(prediction)
}
# Function 3: Inverse Distance Weighting (IDW)
# Weight = 1 / distance^p
# Here p = 2
predict_idw <- function(train, test_point, p = 2) {
distances <- sqrt(
(train$x - test_point$x)^2 +
(train$y - test_point$y)^2
)
# Handle exact coordinate match
if (any(distances == 0)) {
return(train$Moist[which.min(distances)])
}
weights <- 1 / (distances^p)
prediction <- sum(
weights * train$Moist
) / sum(weights)
return(prediction)
}
# Function 4: Ordinary Kriging
# Variogram is estimated using training data only
# to avoid information leakage during cross-validation.
predict_kriging <- function(train, test_point) {
# Convert training data to SpatialPointsDataFrame
coordinates(train) <- ~x+y
# Convert test point to SpatialPoints
test_sp <- data.frame(
x = test_point$x,
y = test_point$y
)
coordinates(test_sp) <- ~x+y
# Initial variogram model
initial_variogram <- vgm(
psill = var(train$Moist),
model = "Exp",
range = max(
1,
max(spDists(
coordinates(train),
longlat = FALSE
)) / 3
),
nugget = 0
)
# Calculate empirical variogram
empirical_variogram <- variogram(
Moist ~ 1,
data = train
)
# Fit variogram model
fitted_variogram <- tryCatch(
fit.variogram(
empirical_variogram,
model = initial_variogram
),
error = function(e) initial_variogram
)
# Ordinary Kriging prediction
kriging_result <- tryCatch(
krige(
Moist ~ 1,
locations = train,
newdata = test_sp,
model = fitted_variogram
),
error = function(e) NULL
)
# Fallback if kriging fails
if (is.null(kriging_result)) {
return(
predict_idw(
as.data.frame(train),
as.data.frame(test_point),
p = 2
)
)
}
return(kriging_result$var1.pred[1])
}
5. LEAVE-ONE-OUT CROSS
VALIDATION
n <- nrow(data_spatial)
# Store predictions
cv_results <- data.frame(
Node_ID = data_spatial$Node_ID,
Observed = data_spatial$Moist,
Closest = NA_real_,
NN = NA_real_,
IDW = NA_real_,
Kriging = NA_real_,
Ensemble = NA_real_
)
for (i in 1:n) {
cat("LOOCV iteration:", i, "of", n, "\n")
# Training data: all points except point i
train_data <- data_spatial[-i, ]
# Test data: point i
test_data <- data_spatial[i, ]
# Closest Observation
pred_closest <- predict_closest(
train_data,
test_data
)
# Nearest Neighbors
pred_nn <- predict_nn(
train_data,
test_data,
k = 3
)
# IDW
pred_idw <- predict_idw(
train_data,
test_data,
p = 2
)
# Ordinary Kriging
pred_kriging <- predict_kriging(
train_data,
test_data
)
# Ensemble of three interpolation methods:
# Closest Observation + IDW + Kriging
pred_ensemble <- mean(
c(
pred_closest,
pred_idw,
pred_kriging
),
na.rm = TRUE
)
# Store predictions
cv_results$Closest[i] <- pred_closest
cv_results$NN[i] <- pred_nn
cv_results$IDW[i] <- pred_idw
cv_results$Kriging[i] <- pred_kriging
cv_results$Ensemble[i] <- pred_ensemble
}
## LOOCV iteration: 1 of 12
## Warning in fit.variogram(empirical_variogram, model = initial_variogram): No
## convergence after 200 iterations: try different initial values?
## [using ordinary kriging]
## LOOCV iteration: 2 of 12
## Warning in fit.variogram(empirical_variogram, model = initial_variogram): No
## convergence after 200 iterations: try different initial values?
## [using ordinary kriging]
## LOOCV iteration: 3 of 12
## Warning in fit.variogram(empirical_variogram, model = initial_variogram): No
## convergence after 200 iterations: try different initial values?
## [using ordinary kriging]
## LOOCV iteration: 4 of 12
## Warning in fit.variogram(empirical_variogram, model = initial_variogram): No
## convergence after 200 iterations: try different initial values?
## [using ordinary kriging]
## LOOCV iteration: 5 of 12
## Warning in fit.variogram(empirical_variogram, model = initial_variogram): No
## convergence after 200 iterations: try different initial values?
## [using ordinary kriging]
## LOOCV iteration: 6 of 12
## Warning in fit.variogram(empirical_variogram, model = initial_variogram): No
## convergence after 200 iterations: try different initial values?
## [using ordinary kriging]
## LOOCV iteration: 7 of 12
## Warning in fit.variogram(empirical_variogram, model = initial_variogram): No
## convergence after 200 iterations: try different initial values?
## [using ordinary kriging]
## LOOCV iteration: 8 of 12
## Warning in fit.variogram(empirical_variogram, model = initial_variogram): No
## convergence after 200 iterations: try different initial values?
## [using ordinary kriging]
## LOOCV iteration: 9 of 12
## Warning in fit.variogram(empirical_variogram, model = initial_variogram): No
## convergence after 200 iterations: try different initial values?
## [using ordinary kriging]
## LOOCV iteration: 10 of 12
## Warning in fit.variogram(empirical_variogram, model = initial_variogram): No
## convergence after 200 iterations: try different initial values?
## [using ordinary kriging]
## LOOCV iteration: 11 of 12
## Warning in fit.variogram(empirical_variogram, model = initial_variogram): No
## convergence after 200 iterations: try different initial values?
## [using ordinary kriging]
## LOOCV iteration: 12 of 12
## Warning in fit.variogram(empirical_variogram, model = initial_variogram): No
## convergence after 200 iterations: try different initial values?
## [using ordinary kriging]
6. CALCULATE RMSE
calculate_rmse <- function(observed, predicted) {
sqrt(
mean(
(observed - predicted)^2,
na.rm = TRUE
)
)
}
rmse_results <- data.frame(
Method = c(
"Closest Observation",
"Nearest Neighbors",
"IDW",
"Kriging",
"Ensemble"
),
RMSE = c(
calculate_rmse(
cv_results$Observed,
cv_results$Closest
),
calculate_rmse(
cv_results$Observed,
cv_results$NN
),
calculate_rmse(
cv_results$Observed,
cv_results$IDW
),
calculate_rmse(
cv_results$Observed,
cv_results$Kriging
),
calculate_rmse(
cv_results$Observed,
cv_results$Ensemble
)
)
)
# Sort from lowest to highest RMSE
rmse_results <- rmse_results %>%
arrange(RMSE)
print(rmse_results)
## Method RMSE
## 1 Nearest Neighbors 8.640270
## 2 Kriging 9.273684
## 3 IDW 9.407563
## 4 Ensemble 10.080028
## 5 Closest Observation 11.974150
7. VISUALIZE RMSE
COMPARISON
ggplot(
rmse_results,
aes(
x = reorder(Method, RMSE),
y = RMSE
)
) +
geom_col() +
coord_flip() +
labs(
title = "Comparison of Interpolation Methods",
subtitle = "Leave-One-Out Cross Validation",
x = "Interpolation Method",
y = "RMSE"
) +
theme_minimal()
8. OBSERVED VS PREDICTED
PLOT
cv_long <- cv_results %>%
pivot_longer(
cols = c(
Closest,
NN,
IDW,
Kriging,
Ensemble
),
names_to = "Method",
values_to = "Predicted"
)
ggplot(
cv_long,
aes(
x = Observed,
y = Predicted
)
) +
geom_point(size = 3) +
geom_abline(
intercept = 0,
slope = 1,
linetype = "dashed"
) +
facet_wrap(~Method) +
labs(
title = "Observed vs Predicted Soil Moisture",
x = "Observed Moisture",
y = "Predicted Moisture"
) +
theme_minimal()
9. SAVE RESULTS
write.csv(
data_spatial,
"~/Downloads/Paraguay_Moisture_Preprocessed.csv",
row.names = FALSE
)
write.csv(
cv_results,
"~/Downloads/Paraguay_LOOCV_Predictions.csv",
row.names = FALSE
)
write.csv(
rmse_results,
"~/Downloads/Paraguay_RMSE_Comparison.csv",
row.names = FALSE
)
cat("Analysis completed successfully!\n")
## Analysis completed successfully!