Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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?

share|improve this question
add comment

2 Answers

up vote 2 down vote accepted

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.

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

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
}
share|improve this answer
    
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? –  Aren Cambre Nov 14 '12 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 '12 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? –  Aren Cambre Nov 14 '12 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 '12 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 '12 at 18:56
show 11 more comments

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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