EsriOceans / btm

Benthic Terrain Modeler
http://resources.arcgis.com/en/communities/oceans/
Mozilla Public License 2.0
26 stars 7 forks source link

Kurtosis and IQR produce invalid results for some edge pixels #128

Closed scw closed 8 years ago

scw commented 8 years ago

For Python 2, one pixel is not being set correctly (position 65, 54 in the raster) and is getting the input bathymetry value for that location. In Python 3, the edge pixels are all retaining their input values, I think this is something related to how the NetCDF layer conversion handles the NetCDF object, see if we can come up with an approach that works along the edges for both environments.

scw commented 8 years ago

Some notes on Kurtosis implementations:

Sample Kurtosis MATLAB documentation

Comparisons:

MATLAB also supports the bias correction (default off)

in_array[0:3,0:3]
# array([[-24.87999916, -24.87000084, -24.95000076],
#       [-24.97999954, -24.88999939, -24.88999939],
#       [-24.96999931, -24.89999962, -24.87999916]], dtype=float32)
# this will by default use Fisher's definition (normal = 0)
scipy.stats.kurtosis(in_array[0:3,0:3].reshape(9,1), axis=0)
# -1.24754453

# now try with Pearson's definition (normal = 3)
scipy.stats.kurtosis(in_array[0:3,0:3].reshape(9,1), fisher=False, axis=0)
# 1.75245547

from noahs' shift grid approach: -1.2475010047310895

converting to float64 first: -1.247544564183044

Now, with R:

library(moments)
window <- (-24.87999916, -24.87000084, -24.95000076, -24.97999954, -24.88999939, -24.88999939, -24.96999931, -24.89999962, -24.87999916)
# "This function computes the estimator of Pearson’s measure of kurtosis."
kurtosis(window)
# 1.752

library(e1071)
window <- (-24.87999916, -24.87000084, -24.95000076, -24.97999954, -24.88999939, -24.88999939, -24.96999931, -24.89999962, -24.87999916)
k
kurtosis(window, type=1)
# -1.248, "type 1", g2 = m4 / m2**2 - 3 == Fisher's
kurtosis(window, type=2)
# -1.233, "type 2", G2 = ((n+1)*g2 + 6) * (n-1)/((n-2)(n-3))
kurtosis(window, type=3)
# -1.615, "type 3", b2 = m4 / s**4 - 3 == (g2+3)(1-1/n)**2 - 3