# Efficient raster sampling of billions of polygons (bounding boxes)

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 :-) – Pimin Konstantin Kefaloukos Mar 5 '12 at 13:13
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? – Pimin Konstantin Kefaloukos Mar 5 '12 at 21:47
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
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?