What is the proper way of calculating the kernel density estimate from geographical coordinates?

Question

I have to calculate the 2d kernel density estimate (kde) from a list of latitude and longitude coordinates. But one degree in latitude is not the same distance as one degree in longitude, this means that the individual kernels would be oval, specially the further the point is from the equator.

In my case the points are all close enough to each other that transforming them to flat earth should not cause many problems. However I'm still curious on how this should be properly handled in case this wasn't true.

Answer 1

You might consider using a kernel especially suitable for the sphere, such as a von Mises-Fisher density

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

where $\mu$ and $\mathbf{x}$ are locations on the unit sphere expressed in 3D Cartesian coordinates.

The analog of the bandwidth is the parameter $\kappa$. The contribution to a location $x$ from an input point at location $\mu$ on the sphere, having weight $\omega(\mu)$, therefore is

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

For each $\mathbf{x}$, sum these contributions over all input points $\mu_i$.

To illustrate, here is R code to compute the von Mises-Fisher density, generate some random locations $\mu_i$ and weights $\omega(\mu_i)$ (12 of them in the code), and display a map of the resulting kernel density for a specified value of $\kappa$ (equal to $6$ in the code).

The points $\mu_i$ are shown as black dots sized to have areas proportional to their weights $\omega(\mu_i)$. The contribution of the large dot near $(100,60)$ is evident throughout the northern latitudes. The bright yellow-white patch around it would be approximately circular when shown in a suitable projection, such as an Orthographic (earth from space).

#
# 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))