You have 2 of the necessary parameters to calculate kurtosis and skewness.
In their simplest form, the Kurtosis is the sum of each value ( x) from the mean ( u) to the 4th power ( x -u)^4, all divided by the variance squared (or standard deviation ^ 4th).. ( sum of ( (x-u)^4) / std dev ^4 ) (skewness (x-u^3, is similarly formed).
Getting the mean and std deviation or variance for each zone is no problem, however, you will still need to get the actual values for each cell in each zone in order to perform the top portion of the equation. This can be accomplished using numpy, but would require separating the values out on a zone basis. Similarly scipy has a kurtosis function but the requirement for single zones would help.
On the upside, you could get some sense of what is going on by getting measures for the other parameter. The sum, range, min, max, majority, minority, and variety etc. could be called into play if you really wanted to try and pull out an estimate of each.
So in short, converting the cells within a zone to points with their associated value could be used. I have some code on the Code Sharing site that does this