11
$\begingroup$

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.

Improve this question
$\endgroup$
6
  • $\begingroup$ As a first guess, I would assume you'd just substitute an appropriate spherical distance metric into a standard kernel approach. $\endgroup$
    – Sycorax ♦
    Sep 16, 2015 at 18:30
  • $\begingroup$ Who's to say that having oval kernels is incorrect? $\endgroup$
    – gung - Reinstate Monica
    Sep 16, 2015 at 18:35
  • 1
    $\begingroup$ @gung Just imagine what would happen if you places a point close enough to a pole. It would be squeezed along the longitudinal axis. And how would you handle a kernel that actually covers one of the poles? $\endgroup$
    – Aaron de Windt
    Sep 16, 2015 at 18:55
  • $\begingroup$ You would have a lump over the pole that is equally high at all longitudes. Why is that incorrect? $\endgroup$
    – gung - Reinstate Monica
    Sep 16, 2015 at 18:56
  • $\begingroup$ @gung Because if I for example choose a kernel diameter of 1 degree it wouldn't be over all longitudes. It would be over 1 longitudinal degree which may be just a few meters if the point is close enough to the pole, compared to the ~110 km that 1 latitudinal degree is. $\endgroup$
    – Aaron de Windt
    Sep 16, 2015 at 19:04

1 Answer 1

Reset to default
7
$\begingroup$

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

[Figure]

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))
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.