adehabitatHR KUD loop not overlaying mapping correctly in R

Question

I have some code which calculates KUD at 95% from fish position data, restricted by the bounds of a shapefile, and then plots this, and overlays the shapefile on top.

########## BOUND KUD TO SHAPEFILE AREA ##########

library(sf)
library(adehabitatHR)
library(raster)
library(viridis)
library(dplyr)

AB_KUD <- read.csv("AB_KUD.csv")

# Get KUD data
coordinates(AB_KUD) <- c("X_UTM","Y_UTM")
proj4string(AB_KUD) <- CRS("+proj=utm +zone=31 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")

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

# Convert KUD to a spatial points dataframe
kud_points <- getverticeshr(kud, percent = 95)

# Convert to raster
r <- raster(extent(Abberton_utm), ncol=100, nrow=100)
r <- rasterize(kud_points, r)

# Read the shapefile
Abberton <- st_read("WFD_Lake_Water_Bodies_Cycle_2.shp")

# Transform to UTM zone 31
Abberton_utm <- st_transform(Abberton, "+proj=utm +zone=31 +datum=WGS84")

# Convert the reservoir boundary to raster
reservoir_raster <- rasterize(Abberton_utm, r)

# Mask the KUD raster using the reservoir boundary raster
kud_masked <- mask(r, reservoir_raster)

# Plot the masked KUD raster
plot(kud_masked, col = viridis(5),
     axes = FALSE, legend = FALSE)  # Adding to the existing plot

axis(1, cex.axis = 2)
axis(2, cex.axis = 2)

# Plot just the geometry of the shapefile
plot(st_geometry(Abberton_utm), col = NA, lwd = 1, border = alpha("black", 0.7), add=TRUE)

I wanted to create a loop which then does this for 90%, then 85% all the way down to 10%. I created the first plot as before, then added a loop that created the KUD data at those values, and plotted it on top of the first plot.

########## BOUND KUD TO SHAPEFILE AREA ##########

library(sf)
library(adehabitatHR)
library(raster)
library(viridis)
library(dplyr)

AB_KUD <- read.csv("AB_KUD.csv")
# Get KUD data
coordinates(AB_KUD) <- c("X_UTM","Y_UTM")
proj4string(AB_KUD) <- CRS("+proj=utm +zone=31 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")

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

# Convert KUD to a spatial points dataframe
kud_points <- getverticeshr(kud, percent = 95)

# Convert to raster
r <- raster(extent(Abberton_utm), ncol=100, nrow=100)
r <- rasterize(kud_points, r)

# Read the shapefile
Abberton <- st_read("WFD_Lake_Water_Bodies_Cycle_2.shp")

# Transform to UTM zone 31
Abberton_utm <- st_transform(Abberton, "+proj=utm +zone=31 +datum=WGS84")

# Convert the reservoir boundary to raster
reservoir_raster <- rasterize(Abberton_utm, r)

# Mask the KUD raster using the reservoir boundary raster
kud_masked <- mask(r, reservoir_raster)

# Plot the masked KUD raster
plot(kud_masked, col = viridis(5),
     axes = FALSE, legend = FALSE)  # Adding to the existing plot

axis(1, cex.axis = 2)
axis(2, cex.axis = 2)

# Add KUD data for different percentiles 
for(x in seq(90,10,-5)){
  # Convert KUD to spatial points dataframe for the current percentile
  kud_points <- getverticeshr(kud, percent = x)
  
  # Convert to raster
  r <- raster(extent(Abberton_utm), ncol=100, nrow=100)
  r <- rasterize(kud_points, r)
  
  # Mask the KUD raster using the reservoir boundary raster
  kud_masked <- mask(r, reservoir_raster)
  
  # Plot the masked KUD raster
  plot(kud_masked, col = viridis(100)[100-x], add = TRUE, legend = FALSE)
}

# Plot just the geometry of the shapefile
plot(st_geometry(Abberton_utm), col = NA, lwd = 1, border = alpha("black", 0.7), add=TRUE)

However for some reason it doesn't overlay perfectly. After the third iteration the plotting shifts so they don't overlay but I don't know why.

I tried to see if the extent or resolution changed for some reason for each iteration but they're the same, so I have no idea what causes this. Any idea idea about why it shifts? All the data to replicate this issue can be found here.

Answer 1

After trying multiple approaches, I could not identify what causes the issue with your code. It is indeed a perplexing situation. I can't say with absolute certainty that this is a result of deprecation, but raster has been creating issues since retirement last year. It may also be an issue with plot(), but I found no definitive information either way.

However, here is a solution that incorporates terra (the successor of raster), and ggplot2. I have, where possible, maintained your current workflow. Some of your object names have been modified to make them explicit e.g. 'kud_points" becomes "sf_poly" as getverticeshr(kud, percent = x) returns a SpatialPolygonsDataFrame, not points. Note also that the 10th percentile area falls outside the mask area, regardless of whether a raster or polygon is used as a mask. I have also switched from using PROJ4 strings to EPSG codes to conform with current convention.

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

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

# Load KUD data, set CRS, and generate SpatialPointsDataFrame 
AB_KUD <- read.csv("AB_KUD.csv")
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)

# Loop percentile vector and add to p
for(x in seq(90, 10, -5)) {
  
  # Convert kud to polygon sf for the current percentile
  sf_poly <- getverticeshr(kud, percent = x) %>%
    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, x)
  
  # Add tmpr to ggplot()
  p <- p +
    geom_spatraster(data = tmpr, show.legend = FALSE)
  
}

# Plot p
p +
  geom_sf(data = Abberton_utm, fill = NA, colour = alpha("black", alpha = 0.7)) +
  scale_fill_gradientn(colours = viridis(100)[seq(95, 10, -5)],
                       na.value = "transparent") +
  coord_sf(datum = 32631) +
  theme(legend.position = "none",
        panel.background = element_blank(),
        panel.border = element_rect(colour = "black", fill = NA, linewidth = 1))