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]},
          {y1: [z21a, z21b, z21c],
           y2: [z22a, z22b, z22c],
           y3: [z23a, z23b, z23c]},
          {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:
  • Are you still trying to find the maximum value of each pixel as in previous questions? – Marc Pfister Apr 8 '14 at 14:38
  • This is a similar script, @MarcPfister. The end product of this script is not to create a raster which contains maximum values. Instead, I need the dictionary to contain the values of the pixels that will be graphed. – Dzinic Apr 8 '14 at 14:40
  • Can you say more about the dictionary? – Marc Pfister Apr 8 '14 at 15:19
  • where is the dictionary, I see only a numpy array – gene Apr 8 '14 at 15:23
  • 1
    If you have the memory, just slice the numpy array instead of building a redundant dictionary. – Marc Pfister Apr 8 '14 at 15:38

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.

| improve this answer | |

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)

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.

| improve this answer | |

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.