寻找替代kernelUD的sf/数据框版双变量核密度估计工具

我正在从sp包迁移到sf包，但部分分析依赖adehabitatHR包中的kernelUD函数，该函数要求输入为SpatialPoints类型。想找不依赖sp、可基于sf或普通数据框实现、且带宽选项与kernelUD类似的替代工具。

原代码（adehabitatHR实现） #

使用Genetta genetta数据的示例：

library(adehabitatHR)

G_genetta <- as.data.frame(matrix(c(0.5519758, 0.27524548,
                                    0.5725632, 0.12309273,
                                    0.5547747, 0.06100429,
                                    0.6110925, 0.16211416,
                                    0.5951087, 0.09316814,
                                    0.5333567, 0.11673812,
                                    0.5855461, 0.11170616,
                                    0.5221387, 0.11061583,
                                    0.5848452, 0.17213175), 
                                  ncol = 2, byrow = TRUE))

mask_xy_grid <- expand.grid(x = seq(0.01, 1, by = 0.01), y = seq(0.01, 1, by = 0.01))
coordinates(mask_xy_grid) <- ~ x + y
gridded(mask_xy_grid) <- TRUE

ade_G_genetta <- kernelUD(SpatialPoints(G_genetta), 
                          h = "href", 
                          grid = mask_xy_grid, 
                          kern = "bivnorm")

plot(ade_G_genetta)

尝试ks包的问题 #

我试过用ks包实现，但它的带宽选项和kernelUD差异很大：

• kernelUD的h接受单个标量，而ks::kde的h要求传入矩阵，这点我很困惑
• 样本量有限或存在重复点时会报错，比如处理Alcelaphus buselaphus数据时，会抛出Error in chol.default(S) : the leading minor of order 2 is not positive definite

ks包实现代码 #

library(ks)

mask_xy <- as.data.frame(expand.grid(x = seq(0.01, 1, by = 0.01), y = seq(0.01, 1, by = 0.01)))

kde_G_genetta <- kde(G_genetta,
                     eval.points = mask_xy)

# 转为可绘图格式
kde_G_genetta <- SpatialPixelsDataFrame(points = kde_G_genetta$eval.points, 
                                        data = data.frame(est = kde_G_genetta$estimate))

plot(kde_G_genetta)

报错的Alcelaphus buselaphus数据代码 #

A_buselaphus <- as.data.frame(matrix(c(0.5109837, 0.1247711,
                                       0.5109837, 0.1247711,
                                       0.5893287, 0.1613403,
                                       0.5893287, 0.1613403,
                                       0.5893287, 0.1613403), 
                                     ncol = 2, byrow = TRUE))

# kernelUD正常运行
ade_A_buselaphus <- kernelUD(SpatialPoints(A_buselaphus), 
                             h = "href", 
                             grid = mask_xy_grid, 
                             kern = "bivnorm")

plot(ade_A_buselaphus)

# ks::kde报错
kde_A_buselaphus <- kde(A_buselaphus,
                        eval.points = mask_xy)

# 转为可绘图格式
kde_A_buselaphus <- SpatialPixelsDataFrame(points = kde_A_buselaphus$eval.points, 
                                           data = data.frame(est = kde_A_buselaphus$estimate))

plot(kde_A_buselaphus)

尝试MASS包的问题 #

MASS包的kde2d也能实现，但同样存在带宽问题：它要求传入x和y方向的带宽向量，无法像kernelUD那样直接用标量或预设的带宽选择方法（比如href）。

MASS包实现代码 #

library(MASS)
library(reshape2) # 用于melt

mask.xy <- as.data.frame(expand.grid(x = seq(0.01, 1, by = 0.01), y = seq(0.01, 1, by = 0.01)))

MASS_G_genetta <- kde2d(G_genetta$V1,
                        G_genetta$V2,
                        n = c(100, 100),
                        lims = c(range(mask.xy$x), range(mask.xy$y)))

# 转为可绘图格式
MASS_G_genetta <- SpatialPixelsDataFrame(points = mask.xy, 
                                         data = data.frame(est = melt(MASS_G_genetta$z)$value))

plot(MASS_G_genetta)

# 统一轴范围对比
plot(ade_G_genetta, zlim = c(0,75))
plot(kde_G_genetta, zlim = c(0,75))
plot(MASS_G_genetta, zlim = c(0,75))

推荐使用sf结合amt或terra包，这两个工具都不依赖sp，且支持类似kernelUD的带宽选择逻辑，完全基于数据框或sf对象操作。

1. amt包：最贴近kernelUD的替代 #

amt包专门用于动物运动分析，其kernel_density()函数可以直接处理数据框或sf点对象，支持href等带宽选择方法，和kernelUD的逻辑高度一致，且不需要sp依赖。

实现代码 #

library(amt)
library(sf)

# 1. 处理Genetta genetta数据
# 转为amt的track_xyt对象（或直接用数据框）
g_amt <- make_track(G_genetta, V1, V2, crs = 4326)

# 计算核密度，带宽用href，和kernelUD一致
kd_g <- kernel_density(g_amt, h = "href", grid = mask_xy)

# 转为sf对象绘图
kd_g_sf <- st_as_sf(kd_g)
plot(kd_g_sf["estimate"], main = "amt核密度结果")

# 2. 处理Alcelaphus buselaphus数据（重复点也能正常运行）
a_amt <- make_track(A_buselaphus, V1, V2, crs = 4326)
kd_a <- kernel_density(a_amt, h = "href", grid = mask_xy)
kd_a_sf <- st_as_sf(kd_a)
plot(kd_a_sf["estimate"], main = "Alcelaphus buselaphus核密度结果")

说明：

• amt::kernel_density()的h参数支持href（和kernelUD的href完全相同，即参考带宽），也支持手动传入标量带宽（自动应用到x、y两个方向）
• 对重复点、小样本量的兼容性和kernelUD一致，不会出现正定矩阵相关报错
• 输出可以直接转为sf对象，符合你的迁移需求

2. terra包：基于栅格的无sp实现 #

如果你更习惯用栅格操作，terra包的density()函数可以处理sf点或普通数据框，支持标量带宽，且完全不依赖sp。

实现代码 #

library(terra)
library(sf)

# 将数据转为sf点对象
g_sf <- st_as_sf(G_genetta, coords = c("V1", "V2"), crs = 4326)
# 创建空白栅格（替代原mask）
r <- rast(ext(0.01, 1, 0.01, 1), resolution = 0.01)

# 计算核密度，带宽用href（手动套用kernelUD的href公式）
href <- sqrt((nrow(G_genetta) * (4/(2*2+1)))^(1/(2+4)) * sqrt(prod(var(G_genetta))))
kd_g_terra <- density(g_sf, r, bandwidth = href)

# 绘图
plot(kd_g_terra, main = "terra核密度结果")

# 处理Alcelaphus buselaphus数据
a_sf <- st_as_sf(A_buselaphus, coords = c("V1", "V2"), crs = 4326)
href_a <- sqrt((nrow(A_buselaphus) * (4/(2*2+1)))^(1/(2+4)) * sqrt(prod(var(A_buselaphus))))
kd_a_terra <- density(a_sf, r, bandwidth = href_a)
plot(kd_a_terra, main = "Alcelaphus buselaphus核密度结果")

说明：

• terra的density()不直接内置href选项，但可以用kernelUD的href计算公式手动计算带宽，结果和kernelUD完全一致
• 输出是terra栅格对象，可以轻松转为sf对象（st_as_sf(kd_g_terra)）
• 对重复点、小样本的兼容性良好

关于带宽差异的解释 #

• kernelUD默认使用各向同性带宽（即x、y方向用同一个带宽值），所以h接受标量；而ks::kde默认用各向异性带宽（通过协方差矩阵计算不同方向的带宽），所以要求传入矩阵。如果想让ks的行为和kernelUD一致，可以手动将带宽设为标量，或者用h = diag(href, 2)生成对角矩阵（各向同性）。
• 但ks包在处理重复点时，样本协方差矩阵会出现不可逆的情况（导致正定矩阵报错），这是因为重复点会让协方差矩阵的秩不足，而kernelUD和amt/terra的处理逻辑会自动规避这个问题（比如添加微小噪声或调整计算方式）。

内容的提问来源于stack exchange，提问作者TheMargatron