从地理坐标计算内核密度估计的正确方法是什么?

机器算法验证 空间的 内核平滑 地质统计学 地理
2022-03-03 01:17:54

我必须从纬度和经度坐标列表中计算 2d 内核密度估计值 (kde)。但是纬度 1 度与经度 1 度的距离不同,这意味着单个内核将是椭圆形的,特别是距离赤道越远的点。

在我的情况下,这些点都彼此足够接近,因此将它们转换为平坦的地球不会引起很多问题。但是,如果这不是真的,我仍然对如何正确处理这件事感到好奇。

1个回答

您可以考虑使用特别适合球体的内核,例如von Mises-Fisher 密度

f(x;κ,μ)exp(κμx)

其中是以 3D 笛卡尔坐标表示的单位球面上的位置。μx

带宽的模拟是参数因此,从球体上处的输入点的贡献,因此是κxμω(μ)

ω(μ)f(x;κ,μ).

对于每个,将这些贡献对所有输入点xμi

为了说明,这里是R计算 von Mises-Fisher 密度的代码,生成一些随机位置和权重(其中 12 个在代码中),并显示指定内核密度的映射的值(在代码中等于)。μiω(μi)κ6

[数字]

显示为大小与其权重成正比的黑点附近的大点的贡献在整个北纬地区都很明显。当以合适的投影(例如正交投影(来自太空的地球))显示时,其周围明亮的黄白色斑块将近似为圆形。μiω(μi)(100,60)

#
# von Mises-Fisher density.
# mu is the location and x the point of evaluation, *each in lon-lat* coordinates.
# Optionally, x is a two-column array.
#
dvonMises <- function(x, mu, kappa, inDegrees=TRUE) {
  lambda <- ifelse(inDegrees, pi/180, 1)
  SphereToCartesian <- function(x) {
    x <- matrix(x, ncol=2)
    t(apply(x, 1, function(y) c(cos(y[2])*c(cos(y[1]), sin(y[1])), sin(y[2]))))
  }
  x <- SphereToCartesian(x * lambda)
  mu <- matrix(SphereToCartesian(mu * lambda), ncol=1)

  c.kappa <- kappa / (2*pi*(exp(kappa) - exp(-kappa)))
  c.kappa * exp(kappa * x %*% mu)
}
#
# Define a grid on which to compute the kernel density estimate.
#
x.coord <- seq(-180, 180, by=2)
y.coord <- seq(-90, 90, by=1)
x <- as.matrix(expand.grid(lon=x.coord, lat=y.coord))
#
# Give the locations.
#
n <- 12
set.seed(17)
mu <- cbind(runif(n, -180, 180), asin(runif(n, -1, 1))*180/pi)
#
# Weight them.
#
weights <- rexp(n)
#
# Compute the kernel density.
#
kappa <- 6
z <- numeric(nrow(x))
for (i in 1:nrow(mu)) {
  z <- z + weights[i] * dvonMises(x, mu[i, ], kappa)
}
z <- matrix(z, nrow=length(x.coord))
#
# Plot the result.
#
image(x.coord, y.coord, z, xlab="Longitude", ylab="Latitude")
points(mu[, 1], mu[, 2], pch=16, cex=sqrt(weights))