Как правильно рассчитать оценку плотности ядра по географическим координатам?

11

Я должен рассчитать оценку плотности ядра (kde) из списка координат широты и долготы. Но один градус по широте - это не то же самое расстояние, что и один градус по долготе, это означает, что отдельные ядра будут овальными, особенно если точка дальше от экватора.

В моем случае все точки достаточно близки друг к другу, поэтому их преобразование в плоскую землю не должно вызывать многих проблем. Однако мне все еще интересно, как это должно быть правильно обработано на случай, если это не так.

Аарон де Виндт
источник
В качестве первого предположения я бы предположил, что вы просто подставили бы подходящую сферическую метрику расстояния в стандартный подход ядра.
Sycorax сообщает, что восстановит Монику
Кто сказал, что овальные ядра - это неправильно?
gung - Восстановить Монику
1
@ gung Только представь, что случится, если ты поставишь точку достаточно близко к полюсу. Было бы сжато вдоль продольной оси. И как бы вы справились с ядром, которое на самом деле покрывает один из полюсов?
Аарон де Виндт
У вас будет шишка над полюсом, которая одинаково высока на всех долготах. Почему это неправильно?
gung - Восстановить Монику
@gung Потому что, если я, например, выберу диаметр ядра 1 градус, это будет не по всем долготам. Было бы более 1 продольного градуса, который может составлять всего несколько метров, если точка находится достаточно близко к полюсу, по сравнению с ~ 110 км, что составляет 1 широтный градус.
Аарон де Винд

Ответы:

7

Вы могли бы рассмотреть использование ядра, особенно подходящего для сферы, такого как плотность фон Мизеса-Фишера

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

где и - местоположения на единичной сфере, выраженные в трехмерных декартовых координатах.μx

Аналогом пропускной способности является параметр . Вклад в местоположение из точки входа в местоположение на сфере, имеющий вес , поэтому равенκxμω(μ)

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

Для каждого суммируйте эти вклады по всем входным точкам .xμi

Чтобы проиллюстрировать это, вот Rкод для вычисления плотности фон Мизеса-Фишера, генерации некоторых случайных местоположений и weights (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))
Whuber
источник