7

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 Apr 10 '18 at 20:00
6

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 Apr 10 '18 at 19:59
  • I am using python 3.6.5 within Anaconda – Kang Apr 10 '18 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. – WhiteboxDev Apr 10 '18 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. – WhiteboxDev Apr 10 '18 at 20:48
3

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)
    else:
        dest = driver.Create(outputpath, width, height, nbands, dtype)

    # Write output raster
    if color_table != 0:
        dest.GetRasterBand(1).SetColorTable(color_table)

    dest.GetRasterBand(1).WriteArray(raster_array)

    if nodata is not False:
        dest.GetRasterBand(1).SetNoDataValue(nodata)

    # Set transform and projection
    dest.SetGeoTransform(gt)
    wkt = data_obj.GetProjection()
    srs = osr.SpatialReference()
    srs.ImportFromWkt(wkt)
    dest.SetProjection(srs.ExportToWkt())

    # 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
    else:
        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 Apr 11 '18 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 Apr 11 '18 at 17:34
  • Thanks very much for your help!!! This is the code I want. – Kang Apr 11 '18 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 Apr 11 '18 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 Jul 23 '18 at 16:41

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.