Calculating zonal entropy index for raster using rasterstats

Question

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:

Answer 1

I actually understood the problem with the rasterstats function right before posting the question. The problem is that it passes an array of arrays to your custom function as input (each array corresponds to a raster column). So, in your custom function you need to flatten the input (I also added a rounding function to keep 2 decimals):

def my_entropy(x):
    return(round(entropy(x.flatten(),base=2),2))

Now each feature has a unique entropy value:

However, I am still interested in solution based on GIS software, if there is any idea.

I also just noticed that the entropy function of scipy must be fed with frequencies, not unique values. So we need to extract the unique values and their frequencies with numpy first, thus my_entropy function becomes:

def my_entropy(x):
    counts=np.unique(x,return_counts=True)[1] #use the return_counts=True and get the array of the counts for each value
    sum_counts=np.sum(counts)
    x=counts/sum_counts
    return(round(entropy(x.flatten(),base=2),2))