Jaki jest właściwy sposób obliczania oszacowania gęstości jądra na podstawie współrzędnych geograficznych?

11

Muszę obliczyć 2d oszacowanie gęstości jądra (kde) z listy współrzędnych szerokości i długości geograficznej. Ale jeden stopień szerokości geograficznej nie jest taki sam, jak jeden stopień długości geograficznej, co oznacza, że ​​poszczególne jądra byłyby owalne, szczególnie im dalej punkt znajduje się od równika.

W moim przypadku wszystkie punkty są wystarczająco blisko siebie, więc przekształcenie ich w płaską ziemię nie powinno powodować wielu problemów. Jednak wciąż jestem ciekawy, jak należy to odpowiednio traktować, na wypadek gdyby to nie była prawda.

Aaron de Windt
źródło
Na pierwszy rzut oka zakładam, że po prostu zastąpisz odpowiednią sferyczną metrykę odległości standardowym podejściem do jądra.
Sycorax mówi Przywróć Monikę
Kto powiedział, że owalne jądra są nieprawidłowe?
gung - Przywróć Monikę
1
@gung Wyobraź sobie, co by się stało, gdybyś umieścił punkt wystarczająco blisko bieguna. Zostałby ściśnięty wzdłuż osi podłużnej. A jak poradziłbyś sobie z jądrem, które faktycznie pokrywa jeden z biegunów?
Aaron de Windt,
Miałbyś guzek nad biegunem, który jest jednakowo wysoki na wszystkich długościach geograficznych. Dlaczego to jest nieprawidłowe?
Gung - Przywróć Monikę
@ gung Bo jeśli na przykład wybiorę średnicę jądra 1 stopień, nie byłoby to na wszystkich długościach geograficznych. Byłoby to ponad 1 stopień podłużny, który może wynosić zaledwie kilka metrów, jeśli punkt jest wystarczająco blisko bieguna, w porównaniu do ~ 110 km, co oznacza 1 stopień podłużny.
Aaron de Windt,

Odpowiedzi:

7

Możesz rozważyć użycie jądra szczególnie odpowiedniego dla kuli, takiego jak gęstość von Misesa-Fishera

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

μx

κxμω(μ)

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

xμi

Rμiω(μi)κ6

[Postać]

μ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
źródło