Python, GDAL, VRT: help wth code efficiency and processing time

Question

I'm going to attempt to make this question as succinct as possible while still including relevant information, as I have trouble with being concise and am trying to avoid writing unnecessary stuff.

Anyway, I need help with increasing the efficiency of my image processing python code. I'm relatively new to python/gdal/coding in general, and don't have very much knowledge about performance, processing time, efficiency, etc. -- so I need your help. However, I do have two specific ideas that could probably speed things up, if I knew how to implement them. In a nutshell, my code does the following:

1. Takes a set of classified Landsat GeoTIFFs (say number of images in set = N), whose sizes all vary by 50-100 pixels, and extracts one element, water, from the classified image. It does this by utilizing GDAL to open and read the classified tiff as an array, and using the where function to produce an array of 1's and 0's (1's where pixel = water and 0's where pixel = not water). Currently, my code writes these arrays as tiffs using the same geospatial information as the original classified tiff, and then uses the set, as a whole, for input of the next step. However, it this step of writing the 'extracted' arrays as tiffs that takes the majority of processing time (8-9 of the 11-12 total minutes per individual N), and so I would like to try a different approach (see bottom) since these 'extracted tiffs' are only needed temporarily for step 2.

2. Next, I use a list pointing to these individual 'extracted' tiffs as the input for gdalbuildvrt, so that a VRT "stack" is built, encompassing the extent of all the individual files. This step is important for implementing step 3 properly, because of the individual file's varying sizes. It is also the main reason why I take the time to write out the 'extracted' arrays as tiffs-- to preserve their spatial information so that gdalbuildvrt will accurately stack each image on top of each other according to their spatial reference. [NOTE: When calling gdalbuildvrt, I use a mask so that the "empty" edges (as in, the empty space around the smaller tiffs) of the individual tiffs are filled with 0's]

3. After creating the N-layer virtual stack of the individual images, and thus 'normalizing' their sizes accurately according to their spatial info, I use a for loop (done since GDAL can only read one band at a time) to read out each resized 'extracted tiff' as a layer of a 3-d numpy array [this array would be of shape (N files, Nlines, Nsamples) if the overall size of the stack is Nlines by Nsamples]. Then, using numpy.sum, I get a 2-d array of size (Nlines, Nsamples), that represents the per-pixel total number of water observations seen in the N files over the course of the set. It then writes this total array as a tiff, which is the final output that I am interested in.

Now knowing my overall goal, hopefully you understand my code explanation. I will attach it if it will be helpful, but it's fairly long and has a lot of comments so I'm leaving it out for now.

There are two ways I have thought about that could (or could not) help improve the efficiency of my code:

1. Skip the step of writing out the 'extracted' arrays as tiffs and either feed them directly into gdalbuildvrt after extracting (while still retaining the spatial information, which it absolutely has to do) or write each file as some other virtual/in-memory format that won't take as long as writing tiffs, which can then be fed as inputs into gdalbuildvrt. I've read that you can't write an array as a VRT, so if that's true, are there any suggestions for alternative methods?

2. I was also thinking of rearranging the code to do this: Feed the classified tiffs directly into gdalbuildvrt which will resize them like it already does, read each classified band one by one into the large 3-d numpy array (in the loop I described above), THEN using the where function to extract water and not water. After that I would total the numpy array, and then write that result to tiff, and be done. The only issue I'm afraid of with this method is running into memory problems when trying to do the where function over a large 3-d numpy array. I imagine the volume of this 3-d numpy array might be the same volume as the numpy arrays in the original method, but I don't know how having values ranging 1-7 (each class) instead of 0's and 1's will affect performance. I also don't know how numpy's where function will act in terms of performance over the large 3-d array.

Can anyone shed some light on what might be the best improvement, in terms of performance and/or efficiency? I realize this is a long post, but I figure being more explicit is better than not enough detail.

Notes: I'm using Python and GDAL for the majority of my code, and cannot resort to arcpy or any other module that requires a license. Also, in case you are wondering, I am going to be processing thousands of 'sets' (~4050) with N ranging anywhere from 10 to 30, so I really need to find ways to improve efficiency. Right now, 11-12 minutes per n puts processing time at a few years... so, I need your help!

Answer 1

I suggest that you build you vrt at the first stage, so that you only process what you need. See this post about lazy computation.

You could also get some inspiration from this post aboutPython function for pixels.

But my last recommendation is to use Orfeo Toolbox (an open source software). You can do both steps with the Band math filter application (a condition, then a sum).

Answer 2

Without seeing your actual script, your best bet is to utilize the full power of your machine. There are several packages that allow systems with multiple processors or cores to run operations in parallel. Some that come to mind are multiprocessing and pp.

If you are associated with a University, it may be worth contacting the computer facility to inquire if you could sign up for time on their super computer or cluster computer. Imagine running 1,000+ operations in parallel!