1

I have a grid of cells (100km x 100km) around the African continent (shapefile) . I am doing a zonal_stats with the population of each cell. It is a raster data (geotiff) with resolution 1kmx1km.

This is the code I am using:

import geopandas as gpd
from pyproj import CRS
import rasterio
from rasterstats import zonal_stats
from rasterio.mask import mask

pop_array, pop_transform = mask(pop, shapes=grid.geometry, crop=True)

zs1 = zonal_stats(grid, pop_array[0], affine=pop_transform, stats=['mean'], nodata=np.nan)
grid['pop'] = [x['mean'] for x in zs1]

where pop is my population raster file and the grid is my shapefile of cells. Both of them have the same crs. However, the following warning appears:

RuntimeWarning: overflow encountered in reduce return umr_sum(a, axis, dtype, out, keepdims, initial, where)

Some of the values of the resulting column "pop" has a value "-inf", and others have a very high negative value, which is surprising for me given that the raster is population count. I have looked at a possible problem and it could be that I need to specify a greater dtype. For example, with astype.

Given that, I have tried to add the following in the last line of my code:

grid['pop'] = [x['mean'] for x in zs1].astype(float128)

But it's still not working with the following error: "AttributeError: 'list' object has no attribute 'astype'"

Any idea about how can I solve it?

"Pop" is the name I assigned to the new column with the results from the zonal stats. It has dtype "float64". When I do pop_array.dtype(), it says to me the following:

TypeError: 'numpy.dtype[float32]' object is not callable

6
  • 1
    What are the dtypes of pop and pop_array?
    – mikewatt
    Feb 28, 2023 at 0:03
  • Thanks @mikewatt for your comment. "Pop" is the name I assigned to the new column with the results from the zonal stats. It has dtype "float64". When I do pop_array.dtype(), it says to me the following: "TypeError: 'numpy.dtype[float32]' object is not callable". I hope this helps.
    – OgeiD
    Mar 1, 2023 at 14:11
  • 1
    I meant the pop var which is being fed into mask(), but I guess that's actually a rasterio dataset. Is the nodata value of that population raster properly defined? If not, you may need to set the nodata arg when calling mask(), so the nodata values get masked out and rasterstats doesn't try to include them when calculating the mean
    – mikewatt
    Mar 1, 2023 at 18:06
  • It works! Also the negative values dissapears. I added "nodata=np.nan" in the mask. If you write the suggestion as an answer, I can give you the bounty. What it is still not clear for me, according to my search in google, the error I had is producing because the number if too big to support the dtype that I have. How is it possible that removing the nodata values, the error dissapear?
    – OgeiD
    Mar 2, 2023 at 8:33
  • And yes, when I do the mask, my pop is already a rasterio dataset, maybe I should have specified that. The raster covers all the world, and I try to cut it just for the size of Africa (grid shapefile)
    – OgeiD
    Mar 2, 2023 at 8:34

1 Answer 1

1
+50

The source nodata value (nan in this case) needs to be defined or supplied as an arg to rasterio.mask.mask(), so in addition to masking out the geometry it also masks out the nodata value and doesn't attempt to include nodata in the calculations.


However I wouldn't expect it to ever overflow, so that's strange. If nans were being included in the calculation the result should be nan as well.

It's possible that somewhere along the line your data was being cast to an integer type. nan is essentially a special-case floating point value which it doesn't apply to integers, so naively treating nan as an int will give you wacky values. It seems like np.mean() is trying to sum all those wacky values and overflowing.

I'm not sure if this is a rasterstats bug or something funky with your specific dataset, like the raster dtype being incorrect or something. To properly debug that I'd probably need your actual input data and the exact versions of all the libraries you're using.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.