Quelle est la bonne façon de calculer l'estimation de la densité du noyau à partir des coordonnées géographiques?

11

Je dois calculer l'estimation de densité de noyau 2d (kde) à partir d'une liste de coordonnées de latitude et de longitude. Mais un degré de latitude n'est pas la même distance qu'un degré de longitude, cela signifie que les grains individuels seraient ovales, spécialement plus le point est éloigné de l'équateur.

Dans mon cas, les points sont tous suffisamment proches les uns des autres pour que leur transformation en terre plate ne pose pas beaucoup de problèmes. Cependant, je suis toujours curieux de savoir comment cela devrait être correctement géré au cas où ce ne serait pas vrai.

Aaron de Windt
la source
À première vue, je suppose que vous substitueriez simplement une métrique de distance sphérique appropriée à une approche de noyau standard.
Sycorax dit Réintégrer Monica le
Qui peut dire qu'avoir des grains ovales est incorrect?
gung - Rétablir Monica
1
@gung Imaginez ce qui se passerait si vous placez un point assez près d'un poteau. Il serait comprimé le long de l'axe longitudinal. Et comment géreriez-vous un noyau qui couvre en fait l'un des pôles?
Aaron de Windt
Vous auriez une bosse sur le poteau qui est également haute à toutes les longitudes. Pourquoi est-ce incorrect?
gung - Rétablir Monica
@gung Parce que si je choisis par exemple un diamètre de noyau de 1 degré, ce ne serait pas sur toutes les longitudes. Il serait supérieur à 1 degré longitudinal, ce qui peut n'être que de quelques mètres si le point est suffisamment proche du pôle, par rapport aux ~ 110 km que représente 1 degré latitudinal.
Aaron de Windt

Réponses:

7

Vous pourriez envisager d'utiliser un noyau particulièrement adapté à la sphère, comme une densité de von Mises-Fisher

F(X;κ,μ)exp(κμX)

μX

κXμω(μ)

ω(μ)F(X;κ,μ).

Xμje

Rμjeω(μje)κ6

[Figure]

μjeω(μje)(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
la source