wafdy_krige.R

library(tidyverse) # wrangling tabular data and plotting
library(sf) # processing spatial vector data
library(sp) # another vector data package necessary for continuity
library(raster) # processing spatial raster data. !!!overwrites dplyr::select!!!
# Load Interpolation Libraries
library(gstat)  # inverse distance weighted, Kriging

da1 <- read.csv("INT_4_ALOS_ASWAN_ZONE_1.csv")
names(da1)[4] <- "dh"

da2 <- read.csv("INT_25_ALOS_ASWAN_ZONE_2.csv")
names(da2)[4] <- "dh"

# Building Grid point to Interpolate Over it
xmin = min(da1$x)
ymin = min(da1$y)
xmax = max(da1$x)
ymax = max(da1$y)
bbox <- c(
  "xmin" = min(da1$x),
  "ymin" = min(da1$y),
  "xmax" = max(da1$x),
  "ymax" = max(da1$y))
grid <- expand.grid(
  x = seq(from = bbox["xmin"], to = bbox["xmax"], by = 30),
  y = seq(from = bbox["ymin"], to = bbox["ymax"], by = 30)
  ) %>% 
  mutate(dh = 0)

# checking the input
all_input <- rbind(
  da1[,c("x", "y", "dh")] %>% mutate(source = "da1"), 
  da2[,c("x", "y", "dh")] %>% mutate(source = "da2")
) %>% 
  rbind(grid %>% mutate(source = "grid"))

# visualize
ggplot(all_input, aes(x, y, color = source, size = source)) + 
  geom_point(alpha = 0.5) +
  scale_size_manual(values = c(5, 5, 0.5) )

##### !!! data of da 2 does not overlap with da1 nor grid!!! #####
# I assume you are trying to krige on da2. We will need another grid for that!

# quick check of input data
ggplot(da2, aes(x, y, color = dh)) + geom_point( size = 3) +
  scale_color_viridis_c()

da2_grid <- expand.grid(
  x = seq(from = min(da2$x), to = max(da2$x), by = 30),
  y = seq(from = min(da2$y), to = max(da2$y), by = 30)
) %>% 
  mutate(dh = 0) %>% 
  st_as_sf(coords = c("x", "y"), crs = 32636, remove = FALSE)

# Dataset da1 does contain only 4 points. I assume you were attempting to krige
# Dataset da2 instead!
ok <- st_as_sf(da2, coords = c("x", "y"), crs = 32636)  

# Automatized Kriging  

ok_result <- krige(
  dh~1, ok, da2_grid,
  model = automap::autofitVariogram(dh~1, as(ok, "Spatial"))$var_model
)

ok_raster <- rasterFromXYZ(
  cbind(
    da2_grid %>% st_drop_geometry() %>% .[,c("x", "y")],
    ok_result$var1.pred
  )
)

plot(ok_raster)