Thursday, 16 February 2017

How to get GDAL to create statistics for GTiff in Python



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:





0
10
5.0189833333333

2.9131294111984




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?



Answer



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:





stats[0]
stats[1]
stats[2]
stats[3]





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


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...