I'd like to re-visit this answer given by Aaron that uses gdal and numpy to calculate the mean value of each band.

Issue 1: With large TIFFs ~ 5-15GB, I am running into a memory error. I am wondering if gdal is uncompressing the TIFF and making it explode?

Issue 2: I switched to a lower resolution and it does run, but the results are different comparing to what QGIS gives in the layer properties. I tired calculating with and without the zero's, but still no match. With zero's is closer, but is still off by about 3 units.

I understand that larger TIFFs will take a bit longer to processes, but it seems faster to plop them into QGIS and get the values from there. Trying to find a "just as fast" automated solution.

I am open to using rasterio, Tiffile, or other raster libraries.

def rgb(ortho):
    raster = gdal.Open(ortho)
    bands = raster.RasterCount

    for band in range(1, bands + 1):
        data = raster.GetRasterBand(band).ReadAsArray().astype('int')
        mean = np.mean(data[data != 0])  # calculate mean without value 0
        print("Band %s: Mean = %s" % (band, round(mean, 2)))


numpy.core._exceptions.MemoryError: Unable to allocate 42.1 GiB for an array with shape (90880, 62208) and data type float64

QGIS Mean Properties:

enter image description here

Python Output with zero's:

enter image description here

Python Output without zero's

enter image description here

  • 2
    Try gdalinfo -stats? Either parsing the output from the command line util, or maybe you can hit that via the python bindings
    – mikewatt
    Nov 2, 2021 at 19:24
  • 3
    With Python use GetStatistics. A full example in thegeoict.com/blog/2020/08/12/….
    – user30184
    Nov 2, 2021 at 19:40
  • 2
    If you want to get accurate statistics you must read every pixel. Gdalinfo supports faster approximated statistics but I do not know how to use it with Python. Perhaps with gdal.Info and InfoOptions gdal.org/python/index.html if you want to try.
    – user30184
    Nov 2, 2021 at 20:31
  • 3
    There seems to be STATISTICS_APPROXIMATE=YES in the QGIS screenshot that you attached. That may explain both the good speed and the difference to the results that you obtained with the numpy method that reads all pixels.
    – user30184
    Nov 2, 2021 at 21:36
  • 4
    Isn't the whole point to get exact statistics? By using the approximated statistics flag you're getting to the wrong answer faster. 15GB isn't particularly large and if they're compressed you can bet your CPU is at max for 1-2 cores but your HDD is barely being tickled, try multi-threading with one process for each core using subprocess popen redirecting stdout to a unique text file per process, it will still take a long time but you'll have many results over the same interval. Are you sure the type is int64, ortho is usually byte, that would help your memory problem. Nov 2, 2021 at 23:31

1 Answer 1


From all the helpful comments and some additional troubleshooting, here is what I found:

QGIS is treating the NaN values (pixels that fall within the extent of the ortho, but are not actually in the ortho - see image) as zeros. So my band averages are darker (closer to zero) then they should be. I sectioned my ortho into multiple pieces and disregarded the sections that contained the edges (for testing purposes) and compared the values from both the original code I posted and from the link user30184 provided to how QGIS calculated their average values and they were negligible.

I also compared the highest res image to a lower res image and found that for my purposes, the mean values were negligible (max difference of around 0.5 units). So for my purpose, I will stick to using the lower res ortho.

Final code:

def rgb_mean(output)
    rgb_df = pandas.DataFrame(columns=['Red', 'Green', 'Blue'])
    raster = gdal.Open(output)
    bands = raster.RasterCount
    avg = []

    for band in range(1, bands + 1):
        data = raster.GetRasterBand(band).ReadAsArray()
        mean = np.mean(data)

    rgb_df = rgb_df.append({'Red': avg[0], 'Green': avg[1], 'Blue': avg[2]}, ignore_index=True)

I removed the [data != 0] from mean = np.mean(data[data != 0] because I do want to include pixels that may have a zero in their band.

enter image description here

  • 1
    Float (32) is the same memory space as int (32), double (64) uses the same memory as int64, both are overkill from a byte data type. Use the GDAL pixel type constants to determine the number of bytes (data type) to use per pixel so you're not wasting memory or degrading your data. In this case I would read in blocks accumulating the values and counting the number of non-nodata values and then divide. Eventually you're going to have a size of a band in uncompressed bytes that is larger than your available RAM and end up thrashing your pagefile which could explain why it's so slow Nov 4, 2021 at 4:27
  • 1
    QGIS and GDAL would leave the NaN values out from the statistics if the value was labeled as nodata. If that was the case then the QGIS image into would show the nodata value in section "bands" and STATISTICS_VALID_PERCENT would be less than 100. You can fix your image with for example gdal_edit.py and assign the nodata value.
    – user30184
    Nov 4, 2021 at 16:09
  • @MichaelStimson I printed my data type print(gdal.GetDataTypeName(data.DataType)) and found it to be byte. To clarify, I can remove .astype('float')?
    – Binx
    Nov 4, 2021 at 18:07
  • 1
    Yes, it should be safe to omit, if you have any problems use astype('byte') to assert the type rather than convert the type. Nov 8, 2021 at 3:41

Your Answer

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

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