I am new when it comes to GIS and relatively new when it comes to programming.

I was trying the “zonal statistics” recipe from the Python GDAL/OGR Cookbook, but I cannot get a decent result. For the code: click here. I know I can also use rasterstats, but I would like to build on this code for some other calculations.

There are 2 problems I am encountering with this code:

  1. It is not taking into account all the features.
  2. It is only calculating the average (and doesn't get any values for mean, median, etc.)

For my test, I am using a raster image of a DEM and a shapefile containing 3 polygons.

Here is my output:

{0: (25.960679824827782, masked, masked, masked, masked), 1: (25.960679824827782, masked, masked, masked, masked), 2: (25.960679824827782, masked, masked, masked, masked)}

I tried printing the top corners and the offsets of every feature it should take into account during the loop, but the problem is that it doesn't change feature (which it should with the GetFeature() method).

I found a similar question as Zonal computations: computing raster pixel counts within polygons using GDAL versus rasterstats, but there was no real answer there.

I am not using arcMap or arcPy

Can anyone help me out?

  • Just use rasterstats like the link you provided suggests. You can modify it to suit your needs--I just did exactly that. I basically have it return all the pixels within my polygon and compute stats in a separate function I wrote. You can use the "raster_out" flag to return masked arrays, too. Not sure what else you'd need. At the very least you can see and understand the parts to build your own.
    – Jon
    Mar 4, 2018 at 5:54
  • Thanks man, I tried rasterstats like you suggested and it provided what I needed anyway. Mar 6, 2018 at 10:27


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.

Browse other questions tagged or ask your own question.