adehabitatHR KUD 循环未在 R 中正确覆盖映射

Question

我有一些代码，可以根据鱼的位置数据计算 95% 的 KUD，受 shapefile 的边界限制，然后绘制此图，并将 shapefile 覆盖在顶部。

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

我想创建一个循环，然后执行 90%、然后 85% 一直到 10%。我像以前一样创建了第一个图，然后添加了一个循环，以这些值创建 KUD 数据，并将其绘制在第一个图的顶部。

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

但是由于某种原因它并没有完美地重叠。第三次迭代后，绘图发生变化，因此它们不会重叠，但我不知道为什么。

我试图查看每次迭代的范围或分辨率是否因某种原因发生变化，但它们是相同的，所以我不知道是什么原因导致的。关于它为什么会改变有什么想法吗？复制此问题的所有数据都可以在here找到。

Answer 1

尝试多种方法后，我无法确定导致代码问题的原因。这确实是一个令人困惑的情况。我不能绝对肯定地说这是弃用的结果，但 raster 和/或

sp

自去年退休以来一直在制造问题。但是，这里有一个包含

terra

（

raster

的后继者）和

ggplot2

的解决方案。在可能的情况下，我保留了您当前的工作流程。我修改了您的一些对象名称以使其明确，例如“kud_points”变为“sf_poly”，因为

getverticeshr(kud, percent = 95)

返回 SpatialPolygonsDataFrame，而不是点。另请注意，无论使用栅格还是多边形作为掩模，第 10 个百分位区域都落在掩模区域之外。

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