Bathymetry[1]

library(sf)

## Warning: package 'sf' was built under R version 4.5.3

## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE

library(sp)
library(gstat)

## Warning: package 'gstat' was built under R version 4.5.3

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.5.3

library(metR)
library(terra)

## Warning: package 'terra' was built under R version 4.5.3

## terra 1.9.25

# Load the data
data <- read.csv("C:/Users/Aldrich Eva/Downloads/Bathymetric Data of Tikub.csv")
lake_shape <- st_read("C:/Users/Aldrich Eva/Downloads/Tikubshapefile.shp", quiet = TRUE)

# Fix lake geometry
sf_use_s2(FALSE)

## Spherical geometry (s2) switched off

lake_shape <- st_make_valid(lake_shape)
lake_shape <- st_union(lake_shape) |> st_as_sf()

if (is.na(st_crs(lake_shape))) {
  st_crs(lake_shape) <- 4326
}

# Preparing sampling points
data$Depth <- abs(data$Depth)

pts <- st_as_sf(
  data,
  coords = c("Longitude", "Latitude"),
  crs = st_crs(lake_shape)
)

pts <- st_intersection(pts, lake_shape)

## although coordinates are longitude/latitude, st_intersection assumes that they
## are planar

## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries

# 4. Project to UTM Zone 51N
lake_utm <- st_transform(lake_shape, 32651)
pts_utm  <- st_transform(pts, 32651)

pts_sp <- as(pts_utm, "Spatial")
lake_v <- vect(lake_utm)

# 5. Create prediction grid
r <- rast(
  ext(lake_v),
  resolution = 3,
  crs = crs(lake_v)
)

grid <- as.data.frame(crds(r))
coordinates(grid) <- ~x + y
proj4string(grid) <- CRS(crs(lake_v))

# 6. Variogram and ordinary kriging
v <- variogram(Depth ~ 1, pts_sp)

depth_var <- var(data$Depth, na.rm = TRUE)

model <- vgm(
  psill = depth_var * 0.7,
  model = "Sph",
  nugget = depth_var * 0.08,
  range = 60
)

fit <- try(fit.variogram(v, model = model), silent = TRUE)

if (inherits(fit, "try-error") || any(is.na(fit$psill)) || any(is.na(fit$range))) {
  fit <- model
}

kr <- krige(Depth ~ 1, pts_sp, grid, model = fit)

## [using ordinary kriging]

# 7. Convert kriging output to raster and mask to lake
kr_df <- as.data.frame(kr)
names(kr_df)[1:2] <- c("x", "y")

r_pred <- rast(
  kr_df[, c("x", "y", "var1.pred")],
  type = "xyz",
  crs = crs(lake_v)
)

r_mask <- mask(r_pred, lake_v)

plot_df <- as.data.frame(r_mask, xy = TRUE, na.rm = TRUE)
names(plot_df)[3] <- "depth"

plot_df$depth[plot_df$depth < 0] <- 0

obs_max <- max(data$Depth, na.rm = TRUE)
plot_df$depth[plot_df$depth > obs_max] <- obs_max

# Create prediction grid
r <- rast(
  ext(lake_v),
  resolution = 3,
  crs = crs(lake_v)
)

grid <- as.data.frame(crds(r))
coordinates(grid) <- ~x + y
proj4string(grid) <- CRS(crs(lake_v))

# Variogram and ordinary kriging
v <- variogram(Depth ~ 1, pts_sp)

depth_var <- var(data$Depth, na.rm = TRUE)

model <- vgm(
  psill = depth_var * 0.7,
  model = "Sph",
  nugget = depth_var * 0.08,
  range = 60
)

fit <- try(fit.variogram(v, model = model), silent = TRUE)

if (inherits(fit, "try-error") || any(is.na(fit$psill)) || any(is.na(fit$range))) {
  fit <- model
}

kr <- krige(Depth ~ 1, pts_sp, grid, model = fit)

## [using ordinary kriging]

# Convert kriging output to raster and mask to lake
kr_df <- as.data.frame(kr)
names(kr_df)[1:2] <- c("x", "y")

r_pred <- rast(
  kr_df[, c("x", "y", "var1.pred")],
  type = "xyz",
  crs = crs(lake_v)
)

r_mask <- mask(r_pred, lake_v)

plot_df <- as.data.frame(r_mask, xy = TRUE, na.rm = TRUE)
names(plot_df)[3] <- "depth"

plot_df$depth[plot_df$depth < 0] <- 0

obs_max <- max(data$Depth, na.rm = TRUE)
plot_df$depth[plot_df$depth > obs_max] <- obs_max

# Contours
minor_breaks <- seq(0, floor(obs_max), by = 5)
major_breaks <- seq(0, floor(obs_max), by = 10)

# Plot
ggplot() +
  geom_tile(
    data = plot_df,
    aes(x = x, y = y, fill = depth),
    width = res(r_mask)[1],
    height = res(r_mask)[2]
  ) +
  geom_contour(
    data = plot_df,
    aes(x = x, y = y, z = depth),
    breaks = minor_breaks,
    color = "grey75",
    linewidth = 0.15
  ) +
  geom_contour(
    data = plot_df,
    aes(x = x, y = y, z = depth),
    breaks = major_breaks,
    color = "black",
    linewidth = 0.3
  ) +
  geom_text_contour(
    data = plot_df,
    aes(x = x, y = y, z = depth),
    breaks = minor_breaks,
    label.placer = label_placer_flattest(),
    size = 2.5,
    colour = "black",
    stroke = 0.1,
    check_overlap = TRUE
  ) +
  geom_sf(
    data = lake_utm,
    fill = NA,
    color = "black",
    linewidth = 0.5
  ) +
  scale_fill_gradientn(
    colours = c(
      "#d9f0ff",
      "#74c0fc",
      "#1c7ed6",
      "#0b3d91",
      "#002458",
      "#000525",
      "#000000"
    ),
    limits = c(0, obs_max),
    breaks = seq(0, floor(obs_max), by = 10),
    name = "depth (m)"
  ) +
  coord_sf(crs = st_crs(lake_utm), expand = FALSE) +
  theme_bw() +
  labs(x = NULL, y = NULL) +
  theme(
    panel.grid = element_blank(),
    aspect.ratio = 1,
    legend.position = "right"
  )