Memory and Processing Efficient Approach to Read a Timeseries in Python

Question

I am trying to store a timeseries of high resolution raster images containing NDVI values into a dictionary for calculations. I am curious as to which approach may the most efficient when loading multiple large files into memory. I currently have tried two options that are both rather inefficient.

1. I would loop through the folder containing all of the images and use GDAL to read the values of the images into a 3D Numpy array.

   img_stack = np.empty((image_width, image_length, len(image_files)), np.dtype('f'))
   
   for i, fname in enumerate(image_files):
       img = gdal.Open(os.path.join(fn + ("\{0}".format(fname)))).ReadAsArray()                 
       img_stack[:, :, i] = img                                                            
       img = None
   

   After being loaded into the array, a loop would go to every single pixel location and retrieve the value of that location in each image.

2. I would do something similar to the previous approach, but instead of loading all of the images into img_stack first, I would go through the images one at a time to reduce memory usage, but it would result in a longer processing time.

Would there be a more efficient way to do this procedure, using Numpy or GDAL?

The structure of the dictionary is as follows:

z_stack = {x1:
          {y1: [z11a, z11b, z11c],
           y2: [z12a, z12b, z12c],
           y3: [z13a, z13b, z13c]},
       x2:
          {y1: [z21a, z21b, z21c],
           y2: [z22a, z22b, z22c],
           y3: [z23a, z23b, z23c]},
       x3:
          {y1: [z31a, z31b, z31c],
           y2: [z32a, z32b, z32c],
           y3: [z33a, z33b, z33c]}}

and the iteration to populate the dictionary is:

z_stack = {}
for x in width:
    z_stack[x] = {}
    for y in length:
        z_stack[x][y] = []
        for z in height:
            z_stack[x][y].append(z)

Answer 1

There's no ultimate answer to this because there's no getting away from memory vs speed tradeoffs.

One way to get 4x "more" memory would be to convert your NDVI arrays from float to signed Int16. You'd still have 4 decimal places precision, and a lot of the float precision was false anyway since the source imagery would have been (I presume) 8 or 16 bit to begin with. Do your math in integer space and convert to floats at the end.

Another way to get more memory would be to distribute the problem over several machines. If your problem allows you to analyze the images one at a time, it should be distributable.

Also, don't overlook Rasterio, https://github.com/mapbox/rasterio, a newer and better Python GDAL library for reading and writing raster data.

Answer 2

I think your first option, memory allowing, is to read all the images into an Ndarray, tuning your storage type as @sgillies says - especially if you're just making graphs. Feed your grapher slices of the array, not new dictionaries.

At the other end of the spectrum, if you're really memory constrained, is to read individual pixel values for each x and y using gdal's ReadRaster or rasterio window slices. In pseudocode something like:

for x in width:
    for y in height:
        series = []
        for image in images:
            i = open(image)
            series.append(i.ReadWindow(x,y,0,0))
        graph(series)

GDAL won't read the whole file just to get one pixel, but it won't be as efficient. Again, it's a tradeoff.

I suppose instead of storing all the time series data, you could store changes. Read images sequentially and store a tuple of the time index and pixel value when the pixel value changes. Again, another tradeoff between computation and memory. Doing some timeit benchmarking would probably help.