Trouble with kerneloverlaphr() after clipping land

Question

I seem to experience a similar issue to the one raised here https://stat.ethz.ch/pipermail/r-sig-geo/2012-December/016858.html, but which never received an answer...

Background:

I have GPS tracking data for multiple seabird colonies and daily fishing effort from Global Fishing Watch. The aim is to compute the kernel UD of all, and then compute the overlap between (a) colonies w/ other colonies, and (b) individual colonies with fishing effort during the time period of the seabird tracking data. For each colony, I have multiple birds and tracks per bird, but as I am interested in the UD of the colony, I treat all individuals as a single "animal".

When I calculate kerneloverlaphr() using the estUD objects returned from kernelUD() function, the outputs from kerneloverlaphr() type = "HR" make sense. However, after I clip the land from my KDEs the overlap values returned are all "1". This makes no sense and can be proven so by plotting the data... I don't understand why kerneloverlaphr() outputs make no sense after clipping the KDEs.

Sample data can be downloaded here: https://www.dropbox.com/sh/yfjy90fh5y9naq0/AADDsfrQ9dPYSOKUFRRcxf2qa?dl=0

My code:

#------------------------------------------------------------------------------
#'*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?