rgeospatialggmapterra

plotting SpatRaster with ggmap in R

Thanks to this post I have been able to calculated some KUD, masked to the area of a shape file, and then plot these with ggplot and SpatRaster from the terra package

However, I would like to map this on to ggmap. My current code is below.

# Create sf object from shapefile, project to WGS84 UTM zone 31N/EPSG:32631
Abberton_utm <- st_read("WFD_Lake_Water_Bodies_Cycle_2.shp") %>%
  st_as_sf() %>%
  st_transform(32631)

#filter data to remove inaccurate positions
AB_filt <- AB_DAT %>% 
  dplyr::filter(HPE <= 2)

#code for plotting KUD data
AB_KUD <- AB_filt %>% 
  dplyr::arrange(Transmitter, DateTime) %>% 
  #dplyr::filter(week == "24") %>% 
  dplyr::select(Species, X_UTM, Y_UTM)

coordinates(AB_KUD) <- c("X_UTM","Y_UTM")
crs(AB_KUD) <- "EPSG:32631"

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

# Create SpatRaster template for rasterize()
# Only need to do this once, no need to replicate it inside your loop
r <- rast(ext(Abberton_utm), ncol = 100, nrow = 100, crs = "EPSG:32631")

# Create 95th pecentile polygon sf from kud for initial ggplot()
sf_poly <- getverticeshr(kud, percent = 95) %>%
  st_as_sf() %>%
  st_set_crs(32631)

# Convert sf_poly to SpatRaster, crop and mask to Abberton_utm
tmpr <- terra::rasterize(sf_poly, r) %>%
  crop(., Abberton_utm, mask = TRUE)

# Assign percentile value to cells in tmpr
tmpr[] <- ifelse(is.na(tmpr[]), NA, 95)

# Create initial plot
p <- ggplot() +
  geom_spatraster(data = tmpr, show.legend = FALSE)

This simply calculates a 95% KUD and uses ggplot to map it. However I'm looking to plot is over a ggmap make with the code below.

# Your code for obtaining the map
Abberton_eastern <- c(0.875, 51.827)
AbRes <- get_map(location=Abberton_eastern, source="google", maptype='stamen_terrain', zoom=14)

I subsequently use a loop that calculates KUD at 5% intervals and overlays them on top of each other in a plot. However, I can't seem to get this first step working. All the data to replicate my code is here.

Solution

  • A few changes to the original question answer are required to get this to work. Primarily, this is to do with the native CRS Google assigns to its data. If you were to project the Google data, it would appear distotrted, especially the text labels. To account for this, all other data will be projected to match the Google data e.g. WGS84/EPSG:4326.

    I've taken some liberties with your stated parameters to create an arguably more balanced map e.g. all of the Abberton Reservoir is included. If this is not desirable, or if you have trouble adapting anything, comment and I will update the answer.

    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())

    result