Qual é a maneira correta de calcular a estimativa da densidade do kernel a partir de coordenadas geográficas?

Eu tenho que calcular a estimativa de densidade de kernel 2d (kde) a partir de uma lista de coordenadas de latitude e longitude. Mas um grau em latitude não é a mesma distância que um grau em longitude, isso significa que os núcleos individuais seriam ovais, especialmente quanto mais longe o ponto estiver do equador.

No meu caso, todos os pontos estão próximos o suficiente um do outro, e transformá-los em terra plana não deve causar muitos problemas. No entanto, ainda estou curioso sobre como isso deve ser tratado adequadamente, caso isso não seja verdade.

spatial kernel-smoothing geostatistics geography Aaron de Windt
fonte

Como primeiro palpite, eu assumiria que você substituísse uma métrica de distância esférica apropriada por uma abordagem padrão do kernel.

Sycorax diz Restabelecer Monica

Quem pode dizer que ter núcleos ovais está incorreto?

gung - Restabelece Monica

@gung Imagine o que aconteceria se você colocasse um ponto próximo o suficiente de um poste. Seria espremido ao longo do eixo longitudinal. E como você lidaria com um kernel que realmente cobre um dos pólos?

Aaron de Windt

Você teria um pedaço sobre o poste que é igualmente alto em todas as longitudes. Por que isso está incorreto?

gung - Restabelece Monica

@gung Porque, por exemplo, se eu escolher um diâmetro de núcleo de 1 grau, não estaria acima de todas as longitudes. Seria acima de 1 grau longitudinal, que pode estar a apenas alguns metros se o ponto estiver próximo o suficiente do poste, em comparação com os ~ 110 km que 1 grau latitudinal é.

Aaron de Windt

Respostas:

Você pode considerar o uso de um kernel especialmente adequado para a esfera, como uma densidade de von Mises-Fisher

f (x; κ, μ) \propto \exp (κ μ^{'} x)

$f(\mathbf{x};\kappa,\mu) \propto \exp(\kappa \mu^\prime \mathbf{x})$

onde e são locais na esfera unitária expressos em coordenadas cartesianas 3D. $\mu$ $\mathbf{x}$

O analógico da largura de banda é o parâmetro . A contribuição para um local de um ponto de entrada no local na esfera, tendo peso , é portanto $\kappa$ $x$ $\mu$ $\omega(\mu)$

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

$\omega(\mu) f(\mathbf{x};\kappa,\mu).$

Para cada , some essas contribuições em todos os pontos de entrada . $\mathbf{x}$ $\mu_i$

Para ilustrar, aqui está o Rcódigo para calcular a densidade de von Mises-Fisher, gerar alguns locais aleatórios e pesos (12 deles no código) e exibir um mapa da densidade resultante do kernel para um determinado valor de (igual a no código). $\mu_i$ $\omega(\mu_i)$ $\kappa$ $6$

Os pontos são mostrados como pontos pretos dimensionados para ter áreas proporcionais aos seus pesos . A contribuição do ponto grande próximo é evidente nas latitudes do norte. A mancha amarela-branca brilhante ao seu redor seria aproximadamente circular quando mostrada em uma projeção adequada, como uma Ortografia (terra do espaço). $\mu_i$ $\omega(\mu_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
fonte