I have an elevation raster and a grid shapefile:
What I want to do is calculate zonal statistics of the raster based on the grid features. The "statistic" I want is the entropy index (Shannon entropy) and it is not included in the tools of most popular software (ArcGIS Desktop and QGIS). I searched for a Python solution and I found the user-defined statistics of rasterstats library: https://pythonhosted.org/rasterstats/manual.html#user-defined-statistics. I attempted the following:
from rasterstats import zonal_stats
from scipy.stats import entropy
def my_entropy(x):
return(entropy(x,base=2))
result=zonal_stats("...path/to/grid.shp","...path/to/elevation.tif", add_stats={'my_entropy':my_entropy})
The code runs, but the result is odd in terms of dimensions. While the length of the result is the same as the number of features of the grid shapefile, the calculated entropy of each feature is an array of values, and not a single value: