How to get GDAL to create statistics for GTiff in Python

Question

I regularly create my own GeoTIFF rasters with GDAL in Python, e.g.:

from osgeo import gdal
from numpy import random
data = random.uniform(0, 10, (300, 200))
driver = gdal.GetDriverByName('GTiff')
ds = driver.Create('MyRaster.tif', 200, 300)
band = ds.GetRasterBand(1)
band.WriteArray(data)
ds = band = None # save, close

however when the result is viewed with ArcCatalog/ArcGIS, it looks either black or grey, since it has no statistics. This is solved either by right-clicking the raster and choosing "Calculate Statistics..." in ArcCatalog (there are several other ways to do this), or using gdalinfo in a command prompt:

gdalinfo -stats MyRaster.tif

will generate MyRaster.tif.aux.xml, which is used by ArcGIS to properly scale the raster. The PAM (Persistent Auxiliary Metadata) file contains the statistics, most notably the minimum and maximum values:

<PAMDataset>
 <PAMRasterBand band="1">
 <Metadata>
 <MDI key="STATISTICS_MINIMUM">0</MDI>
 <MDI key="STATISTICS_MAXIMUM">10</MDI>
 <MDI key="STATISTICS_MEAN">5.0189833333333</MDI>
 <MDI key="STATISTICS_STDDEV">2.9131294111984</MDI>
 </Metadata>
 </PAMRasterBand>
</PAMDataset>

My question: is there a built-in way of getting GDAL to create a statistics file (other than using the gdalinfo -stats command)? Or do I need to write my own?

• Mike T: Same Q/A as gis.stackexchange.com/q/16494

Answer 1

You can use GetStatistics Method to get the stats.

eg.

stats = ds.GetRasterBand(1).GetStatistics(0,1)

it will return (Min, Max, Mean, StdDev)

so the xml can be read:

<PAMDataset>
 <PAMRasterBand band="1">
 <Metadata>
 <MDI key="STATISTICS_MINIMUM">stats[0]</MDI>
 <MDI key="STATISTICS_MAXIMUM">stats[1]</MDI>
 <MDI key="STATISTICS_MEAN">stats[2]</MDI>
 <MDI key="STATISTICS_STDDEV">stats[3]</MDI>
 </Metadata>
 </PAMRasterBand>
</PAMDataset>

I dont know any pythonic way to create/manipulate xml file.But given the simplistic nature of the accompanying xml it should pretty trival to create one it with file I/O operations

• Mike T: It turns out that band.GetStatistics(0,1) will actually calculate the statistics, and add it to the GeoTIFF metadata in the single file. No other files required. However from testing with Esri products, it only works with ArcGIS 10.0 and up, not ArcGIS 9.3 or before.
• user78826: The function is described on the GDAL Page. Based on that, the two arguments passed to the function are bApproxOK (if TRUE statistics may be computed based on overviews or a subset of all tiles) and bForce (if FALSE statistics will only be returned if it can be done without rescanning the image).
• markusN: FWIW, the page is nowadays at gdal.org/api/gdalrasterband_cpp.html

Answer 2

If the statistics are already calculated and included in the file internally, gdalinfo -stats wont create a additional PAM statistics file(.aux.xml) for using GDAL 2.1.0. But its very easy to implement the .xml for your own. Here are some built-in Python modules explained to do that stuff. For myself i used The ElementTree XML API with the code below:

import xml.etree.cElementTree as ET
stats = file.GetRasterBand(band).GetStatistics(0,1)
pamDataset = ET.Element("PAMDataset")
pamRasterband = ET.SubElement(pamDataset, "PAMRasterBand", band="1")
metadata = ET.SubElement(pamRasterband, "Metadata")
ET.SubElement(metadata, "MDI", key = "STATISTICS_MAXIMUM").text = str(stats[1])
ET.SubElement(metadata, "MDI", key = "STATISTICS_MEAN").text = str(stats[2])
ET.SubElement(metadata, "MDI", key = "STATISTICS_MINIMUM").text = str(stats[0])
ET.SubElement(metadata, "MDI", key = "STATISTICS_STDDEV").text = str(stats[3])
tree = ET.ElementTree(pamDataset)
tree.write(destFilePath + ".aux.xml")

The result looks like:

<PAMDataset>
 <PAMRasterBand band="1">
 <Metadata>
 <MDI key="STATISTICS_MINIMUM">-40.65</MDI>
 <MDI key="STATISTICS_MEAN">10.2929293137</MDI>
 <MDI key="STATISTICS_MAXIMUM">45.050012207</MDI>
 <MDI key="STATISTICS_STDDEV">17.4892321447</MDI>
 </Metadata>
 </PAMRasterBand>
</PAMDataset>