How do I compute v.kernel maps in less than 16 hours?

Question

I am using the GRASS GIS 6.4.2 bundled with QGIS 1.8.0, using a PostGIS dataset exported via v.in.org.qgis from QGIS, to create kernel density maps.

The layer has 738,288 points.

My problem is the kernel density map takes about 16 hours to compute.

Here's parameters reported by v.kernel as it's executing:

STDDEV: 1000.000000
RES: 111.419043 ROWS: 458   COLS: 447

Writing output raster map using smooth parameter=1000.000000.

Normalising factor=6482635.018778.

I am running this on a Windows 7 x64 laptop with 8 GB RAM and an Intel Core i7 Q720 1.6 GHz with 4 physical cores. I notice that it's not multithreaded, only using 1 core.

Is there any way to speed this up without harming the quality of the output?

Answer 1

For the record, v.kernel has been improved in GRASS 7 the other day and takes (I believe) now less than 50% of the time.

Answer 2

A kernel density for this size grid only takes a fraction of a second. Evidently, the problem is that v.kernel is processing every one of your three quarters of a million points with too much precision and detail.

Instead, first create a grid to represent the point data, possibly using a finer resolution to reduce the discretization error in location. (Perhaps 4580 rows by 4470 columns, for instance.) The computation even for this much larger grid (it has 100 times as many cells) should still take only a few seconds. Experiment first with a smaller subset of the points and a coarse grid before undertaking the full calculation.

---

Edit

An R FFT in multiple dimensions (tested only in two, though) illustrates the algorithm. It takes roughly three seconds per million cells to convolve moderate-sized kernels with large grids.

filter <- function(x, kernel, ...) {
  # The kernel is centered at its middle.
  # Returns an array of the same dimensions as `x`.
  convolution <- function(x, y) {
    fft.mult <- function(x, ...) {
      z <- x
      d <- dim(z)
      for (i in length(d):1) {
        z <- apply(z, i, fft, ...)
      }
      z
    }
    reverse <- function(y) { # (Can be handy)
      if (is.null(dim(y))) rev(y) else array(rev(y), dim(y))
    }
    x.hat <- fft.mult(x)
    y.hat <- fft.mult(y)
    fft.mult(x.hat * y.hat, inverse=TRUE) / length(x)
  }
  pad <- function(x, pre, post, z=0) {
    # Pad array x in front with pre[i] and in back with post[i] values of z in dimension i.
    d <- length(dim(x))
    if (d > 1) {
      y <- apply(x, d, function(u) pad(u, pre[-d], post[-d], z))
      n <- prod(dim(y)[-d])
      y <- c(rep.int(z, pre[d] * n), y, rep.int(z, post[d] * n))
      array(y, dim(x) + pre + post)
    }
    else {
      y <- c(rep.int(z, pre), x, rep.int(z, post))
      array(y, length(y))
    }
  }
  padding <- dim(kernel)-1
  d <- nextn(dim(x) + padding, ...) # Optional argument is `factors`
  y <- pad(x, padding*0, padding + d - dim(x))
  k <- pad(kernel, padding*0, dim(y) - dim(kernel))
  z <- Re(convolution(y, k))
  e <- dim(x)
  shift <- floor(padding/2)
  z <- do.call(`[`, as.list(c(quote(z), lapply(1:length(e), function(i) 1:e[i]+shift[i]))))
  dim(z) <- dim(x) # Handles mere vectors
  z
}