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:
- It is not taking into account all the features.
- 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?