I am looking for the simplest approach to take the average across multiple raster images, which contain a lot of nan values. The mean value should ignore any nan. The gdal_calc does not work, because it outputs nodata when any nodata value is encoutered.

Is there any simple approach with GDAL or python (with any package)? Prefer not to use ArcGIS.

  • 1
    Is there a reason you can't just load them as arrays using rasterio or gdal and average across the proper axis with np.nanmean? If they're huge, this might not be a viable solution.
    – Jon
    Commented Apr 10, 2018 at 20:00

3 Answers 3


Are you looking to take the average for each grid cell in the stack, or the overall average? If it is the former, you could use the AverageOverlay tool in the WhiteboxTools library. This can be scripted in Python as follows:

from whitebox_tools import WhiteboxTools

wbt = WhiteboxTools()
wbt.work_dir = "/path/to/data/"

wbt.average_overlay(inputs='file1.tif;file2.tif;file3.tif', output='average.tif')

All of the tools in this library will handle NoData values appropriately assuming that the NoData value is flagged properly in the files. For more information on the AverageOverlay tool, see the WhiteboxTools User Manual. For the tool's source code, see here, and in particular, notice how NoData values are handled at this line in the code. That is, the tool will calculate an average for each grid cell only based on the valid values in the image stack. The WhiteboxTools binary executable can be downloaded from the Geomorphometry and Hydrogeomatics Research Group page. There is a WhiteboxTools plugin for QGIS, although it is still somewhat experimental and will be updated soon.

  • Thanks!! This is elegant and simple. I tried , but got the following error: thread 'main' panicked at 'called Result::unwrap()` on an Err value: Os { code: 2, kind: NotFound, message: "No such file or directory" }', libcore/result.rs:945:5 note: Run with RUST_BACKTRACE=1 for a backtrace.`
    – Kang
    Commented Apr 10, 2018 at 19:59
  • I am using python 3.6.5 within Anaconda
    – Kang
    Commented Apr 10, 2018 at 20:04
  • It sounds like a file specification error (as in your file is not being specified correctly in the script). Send me an email (jlindsay at uoguelph .ca) and I'll see if I can get you up and running with it. Commented Apr 10, 2018 at 20:10
  • I think I just figured out what the issue was and it was a simple fix. If you email me, I can get it working for you. I just need to know what operating system you are using. Commented Apr 10, 2018 at 20:48

Since you indicated that reading the rasters into memory would be OK, here's a simple example of how you can achieve your averaging with gdal and numpy. I assume here that all your tiffs are the same size (rows, cols), they are single-banded, share the same CRS, and you do not need to write the computed image to disk.

import gdal
import numpy as np

def write_raster(raster_array, gt, data_obj, outputpath, dtype=gdal.GDT_UInt16, options=0, color_table=0, nbands=1, nodata=False):

    height, width = raster_array.shape

    # Prepare destination file
    driver = gdal.GetDriverByName("GTiff")
    if options != 0:
        dest = driver.Create(outputpath, width, height, nbands, dtype, options)
        dest = driver.Create(outputpath, width, height, nbands, dtype)

    # Write output raster
    if color_table != 0:


    if nodata is not False:

    # Set transform and projection
    wkt = data_obj.GetProjection()
    srs = osr.SpatialReference()

    # Close output raster dataset 
    dest = None

tifflist = ['something1.tif', 'something2.tif'] # the paths to each of the tiffs you want to use in averaging
for i, tiff in enumerate(tifflist):
    gd_obj = gdal.Open(tiff)
    array = gd_obj.ReadAsArray()
    array = np.expand_dims(array,2)
    if i == 0:
        allarrays = array
        allarrays = np.concatenate((allarrays, array), axis=2)
mean_of_tiffs = np.nanmean(allarrays, axis=2)

outputpath = 'wherever you are saving this guy.tif'
write_raster(mean_of_tiffs, gd_obj.GetGeoTransform(), gd_obj, outputpath)
  • Thank you! The code works very well. But I do need to save it out as a geotiff. Sorry that I didn't specify clearly.
    – Kang
    Commented Apr 11, 2018 at 17:18
  • 1
    @Kang added some code to help you save the output array as a geotiff. It will have the same projection and geotransform as the final tiff in your tifflist. I think you should just have to edit the "outputpath" variable and it should run.
    – Jon
    Commented Apr 11, 2018 at 17:34
  • Thanks very much for your help!!! This is the code I want.
    – Kang
    Commented Apr 11, 2018 at 19:47
  • @Kang Note that "write_raster" defaults to a Uint16 datatype; specify gdal.GDT_Float64 if you want floating-point, or you can choose from any of the gdal.GDT_XXXX datatypes. To see them, type gdal.GDT_ in the Spyder console and press "Tab". [Must import gdal first]
    – Jon
    Commented Apr 11, 2018 at 20:28
  • @Jon I followed your code and when I run it I got this error: ValueError: all the input array dimensions except for the concatenation axis must match exactly. I am trying to average tiff files covering a large area where there might be overlap or no overlap to some of them.
    – wondim
    Commented Jul 23, 2018 at 16:41

It's been a while but maybe this is still relevant. I think what you want is the script that follows. Assuming that all rasters are overlapping and have the same size (same tiles, for example), the code below will calculate mean values for each pixel. If there's a NaN value in one or more rasters, it will be ignored and the average value will be calculated from the remaining pixel values.

import rasterio
import numpy as np

def read_raster(name, as_type=float):
    This is a simple function that reads the first band of the raster
    file "name" using rasterio. It is used in the function below.
    r = rasterio.open(name)
    b = r.read(1)
    a = b.astype(as_type)
    return a

def average_rasters(rasters):

    # Empty list to append input raster values
    all_read_rasters = []

    total_input_rasters = len(rasters)

    for r in range(0, total_input_rasters):
        rr = read_raster(rasters[r])

    array_pixel_values = np.array(all_read_rasters)

    # 3D array, calculate along Z axis ("0")
    mean_array = np.nanmean(array_pixel_values, axis=0)
    # Open the first raster in order to get geoinfo!
    temp = rasterio.open(rasters[0])

    with rasterio.open('mean_raster.tif', 'w', driver='Gtiff',
                       width=temp.width, height=temp.height,
                       count=1, crs=temp.crs, transform=temp.transform,
                       dtype='float64') as mean_raster:
        mean_raster.write(mean_array, 1)

# You can add many rasters in the following "source_rasters" list
source_rasters = ['raster1.tif', 'raster2.tif', 'raster3.tif']

# Run the thing!

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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