How can a raster be computed efficiently (in Python), given a set consisting of billions of bounding boxes (read sequentially from a file), and given that the raster values for each cell should give the number of overlapping bounding boxes?

For a 4000 * 4000 raster

I've timed numpy matrix creation:

$ python -m timeit 'import numpy' 'a = numpy.zeros(shape=(4000,4000))'
10 loops, best of 3: 51.7 msec per loop

Standard python matrix creation:

$ python -m timeit 'a = 4000*[0]' 'for i in range(4000):' ' a[i]=4000*[0]'
10 loops, best of 3: 218 msec per loop

So numpy is faster, but still 50 msec per loop, with one billion iterations, yields running time equal to about a year (0.05msec * 1000000000 / 60 / 60 / 24 / 365 = 1.5 years)

So it's not an option to sample each polygon. What is a typical approach for this problem?

  • I want to solve it on a single computer, so no map/reduce solutions please :-) Mar 5 '12 at 13:13
  • 2
    I do not understand the importance of timing raster creation operations. This process needs to create the underlying raster exactly once. Dominating the execution time will be the matter of incrementing counts within the interiors of the bounding boxes. All you have to do is optimize this inner loop. It can be made to go extremely quickly in a compiled language like C or Fortran.
    – whuber
    Mar 5 '12 at 17:46
  • Creating a zero-raster is my crude approximation on how long it would take to increment counts in a bad case. It's a lower bound on how long the worst case takes, where the polygon is as big as the raster, compiled language or not. The real question is, given a 4000x4000 raster, how fast can the entire raster be incremented in C or Fortran on mid-level laptop, back-of-the-envelope? Mar 5 '12 at 21:47
  • 2
    A BB determines a range of rows indexed by i0..i1 and a range of columns j0..j1. In row-by-row storage, you can increment X(i,j0..j1) very rapidly (it's contiguous storage). That can probably be done at around 3E9 increments/sec and even vectorized if you like for much faster operation. Loop i from i0 through i1: that takes care of a single BB. For each BB you have to convert its boundary coordinates into (i0,i1,j0,j1), but that's not much overhead: it can be done faster than you can read the coordinates.
    – whuber
    Mar 5 '12 at 22:56
  • 1
    There is this interesting blog on the ESRI site that talks about using python and multicore processing, may be of help? blogs.esri.com/esri/arcgis/2011/08/29/multiprocessing
    – Hornbydd
    Jun 25 '12 at 20:11

Your timeit includes the numpy import, which would add some overhead. So why don't you write the code for a subset of the bounding boxes and time that loop, then multiply it up to estimate the total running time?

Solving it on a single computer is by its nature serial, and with a relatively simple operation, you might not get any significant optimization from an already simple algorithm. You could try dividing it up in a sort of manual map-reduce operation (I know you have a "no map-reduce" caveat), and running as many instances as you have cores. Mosaicking/merging n rasters (the reduce step) is a trivially fast operation. This will probably be less painful to code than a multi-threaded solution.

Alternatively (or additionally), you could write a program to combine certain bounding boxes such as overlapping or nested ones - this would require a spatial index. If you don't have one, you may find creating one beneficial, especially if you end up locally parallelizing the main algorithm.

Also, don't dismiss multi-computer parallelization out of hand. If your best estimate is over a year, then you need to add up how much money your time will cost in running the single-computer version, and weigh it against hiring some cloud-compute time. As @whuber says, 1024 GPUs will chomp through the data so quickly, it'll cost you next to nothing, even if you spend a week getting your head round CUDA. If it's your boss prohibiting you from trying it on more than one computer, do the cost analysis and hand him some hard numbers - he'll then weigh up the value of the data against the value of your time.


If I understood correctly, what you want is like rendering your set of billions of bounding boxes to an image. Except that instead of "painting" each polygon over a cell (pixel) you count (or accumulate) them.

You can use (relatively) simple code (in OpenGL, Vulcan, Direct3D) to render the polygons and accumulate the count in the stencil buffer. Be careful so that polygons fall on pixel boundaries exactly, and choose a data type for the stencil buffer so that the count does not overflow. I would expect it to run in a few seconds on a single GPU...

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.