I have a stack of rasters (NDVI values over a year) and I want to generate a single raster showing the maximum difference between the layers.

So the result should be a single raster showing the difference between maximum value and minimum value for each point across the landscape. How can I do this in QGIS or Python?

closed as off-topic by ahmadhanb, tinlyx, xunilk, ArMoraer, nmtoken Mar 7 '18 at 10:48

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "Questions seeking help to debug/write/improve code must include the desired behavior, a specific problem or error and the shortest code necessary to reproduce it in the question itself. Providing a clear problem statement and evidence of a code attempt will help others to help you. See: How to create a Minimal, Complete, and Verifiable example." – ahmadhanb, tinlyx, xunilk, ArMoraer, nmtoken
If this question can be reworded to fit the rules in the help center, please edit the question.


If you have GRASS configured with QGIS you can use r.series to do this.

Run r.series from the processing toolbox.

Select your raster layers.

Choose range, and run.

enter image description here

That gives me a new raster in the layer list which overlays my 5 source rasters. If I query these in QGIS I can see the pixel values in my 5 layers and the range pixel value.

enter image description here

In this case r1 and r5 are the raster with the smallest and biggest values, and Aggregate is exactly the difference between them - 5.13 - 1.33


You can use python gdal bindings and numpy to do this.

Here I present a general example:

from osgeo.gdalconst import *
from osgeo import gdal
import numpy as np
import os

arrs = []

for fn in os.listdir(r"folder with rasters"):
   ds = gdal.Open(r"folder with rasters/{}".format(fn))
   arr = np.array(ds.GetRasterBand(1).ReadAsArray())
   del ds

stack = np.stack(arrs)
new_arr = stack.max(axis=0) - stack.min(axis=0)

driver = ds.GetDriver()
dsOut = driver.Create('diff_raster.tif', new_arr.shape[1], new_arr.shape[0], 1, GDT_Float64)
outBand = dsOut.GetRasterBand(1)
outBand.WriteArray(new_arr, 0, 0)



del ds, dsOut, outBand

You can change the new raster Data Type depending on what you need. Just look for an appropiate gdalconst here: http://gdal.org/python/osgeo.gdalconst-module.html/#Variables

  • I get an error when I try the code: Traceback (most recent call last): File "/home/user1/Documents/GIS DATA/SENTINEL2/2017_NDVI_UNCLIPPED/PythonStack/stack.py", line 10, in <module> arr = np.array(ds.GetRasterBand(1).ReadAsArray()) AttributeError: 'NoneType' object has no attribute 'GetRasterBand' – SteveResearch Mar 2 '18 at 12:50
  • It is because gdal did not open the raster. Are you sure you put the right path folder? Usually happens when paths are wrong. – Marcelo Villa Mar 2 '18 at 12:52
  • I added for fn in os.listdir('./clipped_source'): print (fn) and it lists the name of the first file OK but then gives the same "AttributeError: 'NoneType' object has no attribute 'GetRasterBand'" – SteveResearch Mar 2 '18 at 15:56
  • @SteveResearch Yeah, my bad. I just updated the code. fn was only the name so gdal could not find that "path". Try now and let me know if it works – Marcelo Villa Mar 2 '18 at 21:52

Not the answer you're looking for? Browse other questions tagged or ask your own question.