rmapsr-sfkernel-densityadehabitathr

adehabitatHR KUD loop not overlaying mapping correctly in R


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)

enter image description here

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.

enter image description here

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.


Solution

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

    result