web analytics
gdalGeomaticsitopen-sourcePythonRemote SensingTutorial

Getting the statistics of raster with GDAL

GDAL is a powerful geospatial processing library that backups many of the geospatial software you are using. The full list can be found here.

In this post, we will look at a very basic example of reading the raster file including all the bands and looping through the bands to get the statistics (Minimum, Maximum, Mean and Standard Deviation). GDAL has the binding on Python and we will be using Python API to do the same.

    1. The first step is to get the path to your raster file. We will call this data_path.
    2. The next step is to read the raster file. Reading is quite easy with GDAL. Just use the Open method of gdal module. First import the gdal module.
      from osgeo import gdal
      
      ds = gdal.Open(data_path)
      
    3. Next, let’s get the total count of the bands in the raster file. For this, use RasterCount attribute on the ds object obtained above.
      count = ds.RasterCount
      
    4. The raster bands are indexed from 1 to whatever the number of bands are available. To get the raster band (as object and use various methods available to it), we can use GetRasterBand method on the ds object. The argument to the GetRasterBand is the index number of the band. We can then combine the GetStatistics method to the result of GetRasterBand. The first two arguments to the GetStatistics method are the int (or equivalent boolean) values respectively as
      bApproxOK – If TRUE statistics may be computed based on overviews or a subset of all tiles.
      bForce – If FALSE statistics will only be returned if it can be done without rescanning the image.
      For more info on parameteres to the function, see here. The output is the list with Minimum, Maximum, Mean and Standard Deviation resp. for that band in the raster.
      stats = ds.GetRasterBand(1).GetStatistics(0, 1)
      

    Putting it all together

    ds = gdal.Open(data_path)
    count = ds.RasterCount
    response = {}
    for counter in range(1, count+1):
        stats = ds.GetRasterBand(counter).GetStatistics(0, 1)
        response["band_{}".format(counter)] = "Minimum=%.3f, Maximum=%.3f, Mean=%.3f, StdDev=%.3f" % (stats[0], stats[1], stats[2], stats[3])
    print(response)
    

Tags
Show More

Related Articles

Back to top button
Close