剪辑土地后的kerneloverlaphr（）麻烦

Question

我似乎经历了一个类似的问题，在这里提出的问题，但从未收到答案...

背景：

我有多个Seabird colonies colonies colonies colonies colonies的GPS跟踪数据以及全球钓鱼手表的日常捕鱼工作。目的是计算所有人的内核UD，然后计算（a）带有其他菌落的菌落之间的重叠，以及（b）在海鸟跟踪数据的时期内用捕鱼工作的单个菌落。对于每个殖民地，我每只鸟都有多只鸟和轨道，但是当我对殖民地的UD感兴趣时，我将所有人视为一个单一的“动物”。

当我使用从kernelud（）函数返回的estud对象计算kerneloverlaphr（）时，kerneloverlaphr（）type =“ hr”的输出就很有意义。但是，在我从kdes夹住土地后，返回的重叠值都是“ 1”。这是没有意义的，可以通过绘制数据来证明……我不明白为什么剪切KDE后kerneloverlaphr（）输出没有意义。

示例数据可以在此处下载：

我的代码：

#------------------------------------------------------------------------------
#'*Load packages*
#------------------------------------------------------------------------------
library(raster)
library(rworldmap)
library(sp)
library(adehabitatHR)

#-------------------------------------------------------------------------------
#'*1-Create habitat grid for our study area*
#-------------------------------------------------------------------------------

# Make raster layer of study area in LatLong
ras = raster(ext=extent(-70, -50, -60, -38), res=c(0.01,0.01)) #lat/long xmin, xmax, ymin, ymax # global fishing watch data
#give all ponints a "1"
ras[] <- 1
#project the grid to utm or lat/long
projection(ras) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

# load land
worldMap <- getMap(resolution = "high")
projection(worldMap) <- CRS(proj4string(ras))
plot(ras)
plot(worldMap, add = T)

#crop and mask raster by land
#r2 <- crop(ras, worldMap)
r3 <- mask(ras, worldMap, inverse = T, updatevalue=0) ## 1 at sea, 0 on land
plot(r3)
r3.sp <- as(r3, "SpatialPixelsDataFrame")
sea_hab <- r3.sp
image(sea_hab)
proj4string(sea_hab)

#-------------------------------------------------------------------------------
#'*2-load bird tracking data*
#-------------------------------------------------------------------------------

# https://www.dropbox.com/sh/yfjy90fh5y9naq0/AADDsfrQ9dPYSOKUFRRcxf2qa?dl=0 

# Sample data is available here: 
# load seabird data, first use data from Colony A; later use data from Colony B.
df <- read.csv("ColonyA.csv") # use first Colony A, then change to Colony B 
# (This could probably be set up as loop, but I don't know how to make loops - be my guest to suggest one :-D).

# create a spaitialPoints class object with no id.
track_sp <- df[,c("x","y")]
coordinates(track_sp) = ~ x+y # this creates a SpatialPointsDataFrame which will create a KDE for each id.
crs <-"+proj=utm +zone=20 +datum=WGS84 +units=m +no_defs"
proj4string(track_sp) <- crs
track_sp <- spTransform(track_sp, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))

plot(sea_hab)
points(track_sp, pch=".")

#-------------------------------------------------------------------------------
#'*3-Run KDE*
#-------------------------------------------------------------------------------

#  using a kernel density estimate to estimate home range and core range:
track.kde <- kernelUD(track_sp, h = 0.2, grid = sea_hab)

#-------------------------------------------------------------------------------
# 4. now give track_kde a unique and and rerun steps 1 and 2 for next colony
#-------------------------------------------------------------------------------

# Round one: give track.kde a unique name
ColonyA <- track.kde

# Round two
ColonyB <- track.kde

#-------------------------------------------------------------------------------
# 5. now we have two estUD objects; let's compute the overlap 
#-------------------------------------------------------------------------------

# First we need to convert out individual estUDs to a fake estUDm
both_uds <- list(ColonyA=ColonyA, ColonyB=ColonyB)
class(both_uds) <- "estUDm"

# Now we can compute the overlap
kud_overlap <- kerneloverlaphr(both_uds, method = "HR", percent = 50, 
                               0.2, grid = sea_hab)
kud_overlap

## --> this returns values that make sense. 

##########################################

#---------------------------------------------
# 6. Now let's first clip the coastline before computing overlap. So we pick up from step 3.
# Essentially, we now run the code from step 2,3 and 6 for each colony, before giving the output 
# unique names again and calculating the overlap as done in steps 4 and 5.
#---------------------------------------------

# Convert kernel output to a SpatialPixelsDataFrame:
udspdf <- as(track.kde, "SpatialPixelsDataFrame")

# Now convert to a SpatialGridDataFrame, same as the sea habitat map
fullgrid(udspdf) <- TRUE

# Check that the two maps have the same dimensions:
length(udspdf)
length(sea_hab)

# Now calculate the UDs only for the at sea locations.
# For this, we multiply the UD by the habitat variable (1 habitat/ 0 non habitat) and rescale.
resu <- lapply(1:ncol(udspdf), function(i) {
  udspdf[[i]] * sea_hab[[1]] / sum(udspdf[[i]] * sea_hab[[1]])
})

# Now create a dataframe, in which each column represent the UD for a single individual --> this returns just 1 column because we ignore the fact that there are individual animals further up.
resu <- as.data.frame(resu)
names(resu) <- names(udspdf@data)

# Check that the sum of the values is 1
sum(resu)

# Define it as data slot for the udspdf SpatialGridDataFrame
udspdf@data <- resu

# Convert the new udspdf SpatialGridDataFrame to SpatialPixelsDataFrame
fullgrid(udspdf) <- FALSE

# plot it
plot(udspdf)

## We have created a UD for the colony.

## The UD values were calculated for every cell in the grid, but the values on land
## were subsequently set to 0 using the habitat variable in the habitat map.

## For each cell, the associated UD value represents the probability of finding
## an animal there. The sum of the UD values over the whole sea habitat adds up to 1.

## Now we need to transform again to a class-estUDm object to calculate overlap.
# Transform back to estUD
re <- lapply(1:ncol(udspdf), function(i) {
  so <- new("estUD", udspdf[,i])
  so@h <- list(h=0, meth="specified") # specify dummy h values: they are only required to recreate the estUDm
  so@vol <- FALSE
  return(so)
})

names(re) <- names(udspdf)
class(re) <- "estUDm"
image(re)

# now give re the unique colony name, and make it an estUD-class, rather than an estUDm-class object
ColonyAx <- re[[1]]

# and for the second round:
ColonyBx <- re[[1]]

#--------------------------------------
# 6. Ok, now let's repeat steps 5
#--------------------------------------

# First we need to convert out individual estUDs to a fake estUDm
both_uds <- list(ColonyAx=ColonyAx, ColonyBx=ColonyBx)
class(both_uds) <- "estUDm"

# Now we can compute the overlap
kud_overlap <- kerneloverlaphr(both_uds, method = "HR", percent = 50, 
                               0.2, grid = sea_hab)
kud_overlap

# --> this returns only 1s... Why? Why? Why?