Hi,
I am using ArcMap 10.8 version. I have 10 raster files with different data ranging from 0 to 1 for the same study area. I am interested to study how the data from 10 raster files are distributed in each pixel.
I want to calculate the skewness of the data in each pixels.
How can this be done in ArcMap?
Thank you!
Solved! Go to Solution.
better solution,
get your multilayer raster out to a numpy array (RasterToNumPyArray), then use scipy's skew function
a = np.random.random_sample((10, 4, 3)) # 10 layer/band raster with 4 cols 3 rows
from scipy.stats import skew
sk_a = skew(a, axis=0) # Calculate skewness over axis 0 ie local statistic
a0 = a[:, 0, 0] # upper left corner over 10 layers
sk_a0 = skew(a0) # check the skewness for the upper left over the 10 layers
sk_a
array([[-0.55202903, 0.9026473 , -0.22097936],
[ 0.05148204, -0.6512716 , 0.16846484],
[ 0.93081464, 0.00802005, 0.04495503],
[-0.37445255, 0.18772049, -0.10078703]])
sk_a0 # confirm axis 0 calculation
-0.5520290295549968
FYI
You can also bring the results back in. In this example, I created a 10 "band" raster in NumPy and brought it back into ArcGIS Pro
And you can bring it back to numpy
b = arcpy.RasterToNumPyArray(pth)
b.shape
(10, 10, 10)
(a == b).all(-1) # ---- equality check, to make sure input and outputs are equal
array([[1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1]])
Hi Deepa,
If I understand your inputs and analysis objective correctly, you can do this easily using the following workflow. You don't need to create numpy arrays or use any other python package.
Please make sure to use the correct formula for calculating skewness in the map algebra expression below.
Thanks
Noman
Nawajish Noman, PhD
Sr. Principal Product Engineer, Product Engineering Lead
Esri, Redlands, CA.
Input: 10 single band rasters or 1 multiband raster with 10 bands
Steps (using Spatial Analyst tools and operators)
out_skewness = 3 * (cs_mean – cs_median) / cs_stddev
Local statistics Cell Statistics (Spatial Analyst)—ArcGIS Pro | Documentation
can easily give you mean, median and standard deviation.
How difficult the next step would be would depend on which measure of skewness you wanted
Pearson's second skewness coefficient might suffice
better solution,
get your multilayer raster out to a numpy array (RasterToNumPyArray), then use scipy's skew function
a = np.random.random_sample((10, 4, 3)) # 10 layer/band raster with 4 cols 3 rows
from scipy.stats import skew
sk_a = skew(a, axis=0) # Calculate skewness over axis 0 ie local statistic
a0 = a[:, 0, 0] # upper left corner over 10 layers
sk_a0 = skew(a0) # check the skewness for the upper left over the 10 layers
sk_a
array([[-0.55202903, 0.9026473 , -0.22097936],
[ 0.05148204, -0.6512716 , 0.16846484],
[ 0.93081464, 0.00802005, 0.04495503],
[-0.37445255, 0.18772049, -0.10078703]])
sk_a0 # confirm axis 0 calculation
-0.5520290295549968
FYI
You can also bring the results back in. In this example, I created a 10 "band" raster in NumPy and brought it back into ArcGIS Pro
And you can bring it back to numpy
b = arcpy.RasterToNumPyArray(pth)
b.shape
(10, 10, 10)
(a == b).all(-1) # ---- equality check, to make sure input and outputs are equal
array([[1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1]])
Hi Deepa,
If I understand your inputs and analysis objective correctly, you can do this easily using the following workflow. You don't need to create numpy arrays or use any other python package.
Please make sure to use the correct formula for calculating skewness in the map algebra expression below.
Thanks
Noman
Nawajish Noman, PhD
Sr. Principal Product Engineer, Product Engineering Lead
Esri, Redlands, CA.
Input: 10 single band rasters or 1 multiband raster with 10 bands
Steps (using Spatial Analyst tools and operators)
out_skewness = 3 * (cs_mean – cs_median) / cs_stddev