6

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?

2 Answers 2

5

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
}
16
  • I'm not clear what you're saying. First you say to use my 700K points on a finer grid with 100X as many cells, but then you say I should use "a smaller subset of the points and a coarse grid"? Aren't these conflicting? Nov 14, 2012 at 18:23
  • 1
    In writing "experiment first" I mean that for the final calculation you may desire a finer grid, so that you don't lose too much resolution in the point positions, but before waiting for that computation, please make a trial of the complete workflow with fewer points and a coarser grid. This will enable you to streamline the procedure and benchmark its performance so at the end, when you're doing it for real on the full dataset at full resolution, you will know how long it's going to take and what to expect for output.
    – whuber
    Nov 14, 2012 at 18:38
  • OK, thanks. Are you recommending perhaps a random sample of data, then slowly build up the sample size and grid size? But even then, what am I looking for? The right v.kernel parameters to produce the desired output in an acceptable amount of time? The right balance between # of points and grid size? Nov 14, 2012 at 22:33
  • 1
    If all goes well, you will discover that the calculation takes essentially no time at all. You will want to experiment with grid resolution, to assure you don't lose too much precision in gridding the points. You want to make sure that in gridding the points, the data are correctly represented in raster form. You may wish to experiment with the kernel bandwidth--that's usually wise. But don't randomly sample the data: instead, use a small window into one or more characteristic parts of your study area. When all looks good, don't go in stages: just do the final calculation.
    – whuber
    Nov 15, 2012 at 4:11
  • 1
    Increasing the resolution is not intended to make it faster: it's to ensure you aren't losing precision when rasterizing the points. Normally, a good implementation of a density algorithm uses the FFT, which will scale almost with the number of cells in the grid (there's a logarithmic factor in there too, but it's hardly noticeable). A bad implementation will scale with the square of the number of cells. So you should do some timing studies over a range of resolutions with small grids to see what the scaling is. If it's quadratic, look for different software.
    – whuber
    Nov 16, 2012 at 18:56
2

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.

4
  • Is this the same as the v.kernel improvement you put into r53982 and r53983? Nov 24, 2012 at 21:43
  • 2
    Just want to confirm that r53983 of 6.4.3 converted the almost-a-day execution into 25.5 minutes! Nov 25, 2012 at 3:04
  • 1
    ...which is still a thousand times slower than other implementations, such as the R code I supplied!
    – whuber
    Nov 25, 2012 at 16:43
  • Could be - and the quality is the same?
    – markusN
    Nov 26, 2012 at 12:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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