zonal_stats: width and height must be > 0 error

femdias 0
Viewed 30

I am trying to use the function zonal_stats from rasterstats Python package to get the raster statistics from a .tif file of each shape in a .shp file. I manage to do it in QGIS without any problems, but I have to do the same with more than 200 files, which will take a lot of time, so I'm trying the Python way. Both files and replication code are in my Google Drive.

My script is:

import rasterio
import geopandas as gpd
import numpy as np
from rasterio.plot import show
from rasterstats import zonal_stats
from rasterio.transform import Affine


# Import .tif file
raster = rasterio.open(r'M:\PUBLIC\Felipe Dias\Pesquisa\Interpolação Espacial\Arroz_2019-03.tif')

# Read the raster values
array = raster.read(1)

# Get the affine
affine = raster.transform

# Import shape file
shapefile =  gpd.read_file(r'M:\PUBLIC\Felipe Dias\Pesquisa\Interpolação Espacial\Setores_Censit_SP_WGS84.shp')

# Zonal stats
zs_shapefile = zonal_stats(shapefile, array, affine = affine,
                        stats=['min', 'max', 'mean', 'median', 'majority'])

I get the following error:

Input In [1] in <cell line: 22>
    zs_shapefile = zonal_stats(shapefile, array, affine = affine,

  File ~\Anaconda3\lib\site-packages\rasterstats\main.py:32 in zonal_stats
    return list(gen_zonal_stats(*args, **kwargs))

  File ~\Anaconda3\lib\site-packages\rasterstats\main.py:164 in gen_zonal_stats
    rv_array = rasterize_geom(geom, like=fsrc, all_touched=all_touched)

  File ~\Anaconda3\lib\site-packages\rasterstats\utils.py:41 in rasterize_geom
    rv_array = features.rasterize(

  File ~\Anaconda3\lib\site-packages\rasterio\env.py:387 in wrapper
    return f(*args, **kwds)

  File ~\Anaconda3\lib\site-packages\rasterio\features.py:353 in rasterize
    raise ValueError("width and height must be > 0")

I have found this question about the same problem, but I can't make it work with the solution: I have tried to reverse the signal of the items in the Affine of my raster data, but I couldn't make it work:

''' Trying to use the same solution of question:  https://stackoverflow.com/questions/62010050/from-zonal-stats-i-get-this-error-valueerror-width-and-height-must-be-0    '''

old_tif = rasterio.open(r'M:\PUBLIC\Felipe Dias\Pesquisa\Interpolação Espacial\Arroz_2019-03.tif')
print(old_tif.profile) # copy & paste the output and change signs 

new_tif_profile = old_tif.profile

# Affine(0.004611149999999995, 0.0, -46.828504575,
#        0.0, 0.006521380000000008, -24.01169169)


new_tif_profile['transform'] = Affine(0.004611149999999995, 0.0, -46.828504575,
                                      0.0, -0.006521380000000008, 24.01169169)

new_tif_array = old_tif.read(1)
new_tif_array = np.fliplr(np.flip(new_tif_array))

with rasterio.open(r'M:\PUBLIC\Felipe Dias\Pesquisa\Interpolação Espacial\tentativa.tif', "w", **new_tif_profile) as dest:
    dest.write(new_tif_array, indexes=1)


dem = rasterio.open(r'M:\PUBLIC\Felipe Dias\Pesquisa\Interpolação Espacial\tentativa.tif')

# Read the raster values
array = dem.read(1)

# Get the affine
affine = dem.transform

# Import shape file
shapefile =  gpd.read_file(r'M:\PUBLIC\Felipe Dias\Pesquisa\Interpolação Espacial\Setores_Censit_SP_WGS84.shp')

# Zonal stats
zs_shapefile = zonal_stats(shapefile, array, affine=affine,
                        stats=['min', 'max', 'mean', 'median', 'majority'])

Doing this way, I don't get the "width and height must be > 0" error! But every stat in zs_shapefile is "NoneType", so it doesn't help my problem.

Does anyone understands why this error happens, and which sign I have to reverse for making it work? Thanks in advance!

gdal raster python
1 Answers

I would be careful with overriding the geotransform of your raster like this, unless you are really convinced the original metadata is incorrect. I'm not too familiar with Affine, but it looks like you're setting the latitude now as positive? Placing the raster on the northern hemisphere. My guess would be that this lack of intersection between the vector and raster causes the NoneType results.

I'm also not familiar with raster_stats, but I'm guessing it boils down to GDAL & Numpy at the core of it. So something you could try as a test is to add the all_touched=True keyword:
https://pythonhosted.org/rasterstats/manual.html#rasterization-strategy

If that works, it might indicate that the rasterization fails because your polygons are so small compared to the pixels, that the default rasterization method results in a rasterized polygon of size 0 (in at least one of the dimensions). And that's what the error also hints at (my guess).

Keep in mind that all_touched=True changes the stats you get in result, so I would only do it for testing, or if you're comfortable with this difference.

If you really need a valid value for these (too) small polygons, there are a few workarounds you could try. Something I've done is to simply take the centroid for these polygons, and take the value of the pixel where this centroid falls on.

A potential way to identify these polygons would be to use all_touched with the "count" statistic, every polygon with a count of only 1 might be too small to get rasterized correctly. To really find this out you would probably have to do the rasterization yourself using GDAL, given that raster_stats doesn't seem to allow it.

Note that due to the shape of some of the polygons you use, the centroid might fall outside of the polygon. But given how course your raster data is, relative to the vector, I don't think it would impact the result all that much.

An alternative is, instead of modifying the vector, to significantly increase the resolution of your raster. You could use gdal_translate to output this to a VRT, with some form of resampling, and avoid having to write this data to disk. Once the resolution is high enough that all polygons rasterize to at least a 1x1 array, it should probably work. But your polygons are tiny compared to the pixels, so it'll be a lot. You could guess it, or analyze the envelopes of all polygons. For example take the smallest edge of the envelope as more or less the resolution that's necessary for a correct rasterization.

Edit; To clarify the above a bit further.

The default rasterization strategy of GDAL (all_touched=False) is to consider a pixel "within" the polygon if the centroid of the pixel intersects with the polygon.

Using QGIS you can for example convert the pixels to points, and then do a spatial join with your vector. If you remove polygons that can't be joined (there's a checkbox), you'll get a different vector that most likely should work with raster_stats, given your current raster.

You could perhaps use that in the normal way (all_touched=False), and get the stats for the small polygons using all_touched=True.

In the image below, the green polygons are the ones that intersect with the centroid of a pixel, the red ones don't (and those are probably the ones raster_stats "tries" to rasterize to a size 0 array).

enter image description here

Rutger Kassies0
Related