Wrong mean/std stats

[yurithefury]

Hi,

I think there is an error in calculating mean and standard deviation in rasterstats 0.10.3
Here is the statistics calculated in ArcMap for my region:
Band COUNT MIN MAX MEAN STD MEDIAN
1 151057 1378 2399 1584 71 1577
2 151057 1025 2753 1318 119 1305
3 151057 594 3271 1035 191 1014
4 151057 380 3636 854 223 826
5 151057 113 4225 753 252 718
6 151057 366 4869 1853 642 1808
7 151057 153 7181 2947 1099 2873
8 151057 121 5926 2577 916 2519
9 151057 13880 59141 51615 3675 52475
10 151057 14939 44531 37989 1722 38133
11 151057 10063 56708 48523 3707 49289
12 151057 13024 46146 34272 5463 34959
13 151057 23169 59005 46253 3292 46804

And here is the same statistics calculated using rasterstats (have a look at mean and std values for bands 9-13):
Band count min max mean std percentile_50
1 151057 1378 2399 1584 71 1577
2 151057 1025 2753 1318 119 1305
3 151057 594 3271 1035 191 1014
4 151057 380 3636 854 223 826
5 151057 113 4225 753 252 718
6 151057 366 4869 1853 642 1808
7 151057 153 7181 2947 1099 2873
8 151057 121 5926 2577 916 2519
9 151057 13880 59141 23182 28669 52475
10 151057 14939 44531 9556 28485 38133
11 151057 10063 56708 20090 28673 49289
12 151057 13024 46146 5839 28953 34959
13 151057 23169 59005 17820 28623 46804

The code to reproduce the results using rasterstats:
for i in range(1, 14):
stats = zonal_stats('F:/test.shp', "F:/test_2.dat", stats=['count','min', 'max', 'mean', 'std', 'percentile_50'], all_touched=False, band_num=i, nodata=0)

Here is the data: https://www.dropbox.com/sh/e8g0jwr1ci0vp2r/AAATlwU2eFjtbtMCO3R5SB1Fa?dl=0

Can you please let me know what the problem could be?

[yurithefury]

UPDATE: when I convert the image from 16 bit unsigned to 32 bit float, mean and std are calculated correctly

[perrygeo]

I'm not able to reproduce. Even with the uint16 image (test_2.dat), I get results identical to the ArcMap results with both rasterstats==0.10.3 and current master branch:

Band COUNT MIN MAX MEAN STD MEDIAN
1 151057 1378 2399 1584 71 1577
2 151057 1025 2753 1318 119 1305
3 151057 594 3271 1035 191 1014
4 151057 380 3636 854 223 826
5 151057 113 4225 753 252 718
6 151057 366 4869 1853 642 1808
7 151057 153 7181 2947 1099 2873
8 151057 121 5926 2577 916 2519
9 151057 13880 59141 51615 3675 52475
10 151057 14939 44531 37989 1722 38133
11 151057 10063 56708 48523 3707 49289
12 151057 13024 46146 34272 5463 34959
13 151057 23169 59005 46253 3292 46804

I can't see how your second table was produced, unless it was using a different raster than the one you provided in the dropbox link.

[yurithefury]

Thanks for looking into this @perrygeo
Unfortunately my issue remains, and it's definitely not related to the use of a different raster. To make our comparisons easier I extracted just one band (uint16) (see: http://goo.gl/XYbIcJ) where after running:
stats = zonal_stats('test.shp', 'band1.dat', stats="mean", nodata=0)
I'm getting the mean of 9556 instead of 37989. And as I said before, after converting uint16 to float32 I'm getting the correct result.
I guess, it has something to do with my setup. I'm testing on Windows x64, and I installed rasterstats using rasterstats-0.10.3-py2-none-any.whl. The funny thing is that I also ran the same code under Ubuntu 14.04 using rasterstats==0.6.2 and I'm getting the correct result (i.e. 37989).
I guess I am forced to convert all my imagery into float32 format now (under Windows x64).

[sgoodm]

what version of rasterio are you using? try updating that if needed

[yurithefury]

It's up-to-date (i.e. rasterio==0.36.0) and was installed using rasterio-0.36.0-cp27-cp27m-win_amd64.whl

[perrygeo]

This looks like a an issue specific to your platform/installation. I can't reproduce it on any installations and its not clear that there is anything I can do from the code side. My best guess is some sort of binary incompatibility between some of the C components - improperly installed versions of numpy, GDAL, rasterio, shapely, fiona and python itself (maybe its 32-bit?) might be at fault here.

If you can find a way to reproduce this error in another environment, please reopen.

[jonschwenk]

I am having a very similar issue. My raster data are uint16. QGIS computes what I believe are correct stats, while rasterstats gives garbage values (for sum). Unfortunately, my raster is global extent, 30m so converting it to float isn't feasible (and I have multiple global rasters in uintXX). Also on a Win10 x64 system. Rasterstats version 0.12.0, rasterio version 0.36.0.

[jonschwenk]

I should note that I was able to debug this (not sure if it needs to be "fixed" in more places, though).

When you create the masked array, numpy by default makes it the same datatype as the input data. When I added ", dtype=np.float" to line 170 here, I got the same values as QGIS (and much much faster, thanks for writing this software!). I'm not sure why this behavior exists; must have something to do with the way numpy stores intermediate results when calling sum()?

[perrygeo]

@jeueblah thanks for the report. Given that we have multiple reports of the problem on win10 x64, I'm going to reopen and investigate further.

I don't currently have access to a Win 10 machine but I'll do my best.

[perrygeo]

I think we might have a numpy distribution issue.

❓ Can anyone reproduce the problem using:

• windows 64bit
• Official Python.org distribution
• numpy 1.14 wheels from pypi

[jonschwenk]

I am using numpy 1.14 also, but installed within an anaconda environment from their repository. Also, if it helps, I was getting very large negative numbers (~2^31) as the output for sum.

[perrygeo]

I found a windows machine to test on :-) but I don't have a copy of the data to reproduce it and the dropbox links are 404ing :-(

@Fri-end can you send me an updated link to the dataset?

[yurithefury]

@perrygeo unfortunately that dataset is long gone