Apa cara yang tepat untuk menghitung estimasi kepadatan kernel dari koordinat geografis?

11

Saya harus menghitung estimasi kepadatan kernel 2d (KDE) dari daftar koordinat lintang dan bujur. Tetapi satu derajat dalam garis lintang tidak sama dengan satu derajat dalam garis bujur, ini berarti bahwa masing-masing kernel akan berbentuk oval, khususnya titik selanjutnya dari khatulistiwa.

Dalam kasus saya, semua titik cukup dekat satu sama lain sehingga mengubahnya menjadi tanah datar seharusnya tidak menimbulkan banyak masalah. Namun saya masih penasaran tentang bagaimana ini harus ditangani dengan benar jika ini tidak benar.

Aaron de Windt
sumber
Sebagai tebakan pertama, saya akan menganggap Anda hanya mengganti metrik jarak bola yang sesuai menjadi pendekatan kernel standar.
Sycorax berkata Reinstate Monica
Siapa bilang tidak memiliki kernel oval?
gung - Reinstate Monica
1
@ung Bayangkan saja apa yang akan terjadi jika Anda menempatkan titik yang cukup dekat ke sebuah tiang. Itu akan terjepit di sepanjang sumbu longitudinal. Dan bagaimana Anda menangani kernel yang sebenarnya menutupi salah satu kutub?
Aaron de Windt
Anda akan memiliki benjolan di atas tiang yang sama tingginya di semua garis bujur. Mengapa itu salah?
gung - Reinstate Monica
@ung Karena karena saya misalnya memilih diameter kernel 1 derajat tidak akan lebih dari semua garis bujur. Itu akan menjadi lebih dari 1 derajat longitudinal yang mungkin hanya beberapa meter jika titiknya cukup dekat dengan kutub, dibandingkan dengan ~ 110 km seperti 1 derajat lintang.
Aaron de Windt

Jawaban:

7

Anda mungkin mempertimbangkan untuk menggunakan kernel yang sangat cocok untuk bola, seperti kepadatan von Mises-Fisher

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

di mana dan x adalah lokasi pada unit sphere yang dinyatakan dalam koordinat Cartesian 3D.μx

Analog dari bandwidth adalah parameter . Kontribusi ke lokasi x dari titik input di lokasi μ pada bola, memiliki bobot ω ( μ ) , oleh karena itu adalahκxμω(μ)

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

Untuk setiap , jumlah kontribusi ini pada semua titik input μ i .xμi

Untuk menggambarkan, di sini adalah Rkode untuk menghitung kepadatan von Mises-Fisher, menghasilkan beberapa lokasi acak dan bobot ω ( μ i ) (12 dari mereka dalam kode), dan menampilkan peta kepadatan kernel yang dihasilkan untuk nilai tertentu dari κ (sama dengan 6 dalam kode).μiω(μi)κ6

[Angka]

μ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
sumber