Построение SpatRaster с помощью ggmap в R

Чтобы это заработало, необходимо внести несколько изменений в исходный ответ на вопрос. В первую очередь это связано с собственной CRS, которую Google присваивает своим данным. Если бы вы проецировали данные Google, они выглядели бы искаженными, особенно текстовые метки. Чтобы учесть это, все остальные данные будут совпадать с данными Google, например. WGS84/EPSG: 4326.

Я позволил себе некоторые вольности с указанными вами параметрами, чтобы создать, возможно, более сбалансированную карту, например. включено все водохранилище Аббертон. Если это нежелательно или у вас возникли проблемы с адаптацией, прокомментируйте, и я обновлю ответ.

library(sf)
library(adehabitatHR)
library(raster)
library(viridis)
library(dplyr)
library(terra)
library(tidyterra)
library(ggmap)
library(ggplot2)

# Register Google API key
register_google(key = "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX")

# Set sf EPSG CRS value native to map source e.g. ggmap == 4326
mapcrs <- 4326

# Set SpatRaster CRS
mapcrsr <- "EPSG:4326"

# Create sf object from shapefile
Abberton <- st_read("WFD_Lake_Water_Bodies_Cycle_2.shp") %>%
  st_as_sf() %>%
  st_transform(mapcrs)

# Create sf point object, set data crs, project to mapcrs, convert to sp
AB_KUD <- read.csv("AB_KUD.csv") %>%
  st_as_sf(coords = c("X_UTM","Y_UTM"), crs = 32631) %>%
  st_transform(mapcrs) %>%
  as_Spatial()

# Create estUDm object
kud <- kernelUD(AB_KUD[,1], h = "href")  # href = the reference bandwidth

# Create SpatRaster template for rasterize()
r <- rast(ext(Abberton), ncol = 100, nrow = 100, crs = mapcrsr)

# Function to get polygon from bounding box (basis for plot extent and Google map)
bbox_poly <- function(x) {
  
  bb <- sf::st_bbox(x)
  
  p <- matrix(
    c(bb["xmin"], bb["ymin"], 
      bb["xmin"], bb["ymax"],
      bb["xmax"], bb["ymax"], 
      bb["xmax"], bb["ymin"], 
      bb["xmin"], bb["ymin"]),
    ncol = 2, byrow = T
    
  )
  
  sf::st_polygon(list(p))
  
}

# Create box sf
box <- st_sfc(bbox_poly(Abberton)) %>%
  st_set_crs(mapcrs)

# Return st_bbox() coordinates for plot extent
mapext <- st_bbox(box)

# Get box centroid for get_googlemap() centre
st_coordinates(st_centroid(box))
#              X        Y
# [1,] 0.8564451 51.82594

# Define map centre
Abberton_eastern <- c(lon = 0.8564451, lat = 51.82594)

# Get map from Google API
tmp <- get_googlemap(center = Abberton_eastern, # US spelling only!
                     source = "google",
                     maptype = "terrain",
                     scale = 2,
                     zoom = 13,
                     style = "feature:poi|visibility:off") # Remove all POI pins

# Convert Google map object to SpatRaster                     
AbRes <- rast(tmp)

# Write rast() to working directory
writeRaster(AbRes, "AbRes_z13_scale2.tif", overwrite = TRUE)

# Create initial plot
p <- ggplot() +
  geom_spatraster_rgb(data = AbRes)

# Loop percentile vector and add to p
for(x in seq(95, 10, -5)){
  
  # Convert kud to polygon sf for the current percentile
  sf_poly <- getverticeshr(kud, percent = x) %>%
    st_as_sf() %>%
    st_set_crs(mapcrs)
  
  # Convert sf_poly to SpatRaster, crop and mask to Abberton_utm
  tmpr <- terra::rasterize(sf_poly, r) %>%
    crop(., Abberton, mask = TRUE)
  
  # Assign percentile value to cells in tmpr
  tmpr[] <- ifelse(is.na(tmpr[]), NA, x)
  
  # Add tmpr to ggplot()
  p <- p +
    geom_spatraster(data = tmpr, show.legend = FALSE)
  
}

# Set plot extent offset in map units
extoff <- 0.005

# Plot p
p +
  # geom_sf(data = Abberton, fill = NA, colour = alpha("black", alpha = 0.7)) +
  scale_fill_gradientn(colours = viridis(100)[seq(95, 10, -5)],
                       na.value = "transparent") +
  coord_sf(expand = FALSE,
           xlim = c(mapext[[1]] - extoff, mapext[[3]] + extoff),
           ylim = c(mapext[[2]] - extoff, mapext[[4]] + extoff)) +
  theme(legend.position = "none",
        panel.background = element_blank())