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