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*' 'for i in range(4000):' ' a[i]=4000*' 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?