1

I am using the function zonal_stats from the Rasterstats library. I have already used this function on data I have for precipitation which works impeccably. But when i try to run the function using the same vector, but with a different raster dataset (for Actual Evapotranspiration) I get the error message: 

    Traceback (most recent call last):
  File "<input>", line 3, in <module>
  File "/Users/ida/opt/anaconda3/envs/thesis_env/lib/python3.7/site-packages/rasterstats/main.py", line 31, in zonal_stats
    return list(gen_zonal_stats(*args, **kwargs))
  File "/Users/ida/opt/anaconda3/envs/thesis_env/lib/python3.7/site-packages/rasterstats/main.py", line 159, in gen_zonal_stats
    rv_array = rasterize_geom(geom, like=fsrc, all_touched=all_touched)
  File "/Users/ida/opt/anaconda3/envs/thesis_env/lib/python3.7/site-packages/rasterstats/utils.py", line 47, in rasterize_geom
    all_touched=all_touched)
  File "/Users/ida/opt/anaconda3/envs/thesis_env/lib/python3.7/site-packages/rasterio/env.py", line 386, in wrapper
    return f(*args, **kwds)
  File "/Users/ida/opt/anaconda3/envs/thesis_env/lib/python3.7/site-packages/rasterio/features.py", line 347, in rasterize
    raise ValueError("width and height must be > 0")

The code I try to run is this:

vec = '/path/to/file/watersheds_template.shp'
AET = '/path/to/file/AET_2003.tif'
avg_AET = zonal_stats(vec, AET, layer='watersheds_template', stats='mean', geojson_out=True)

I have opened the AET raster file in Qgis, where there is no problem, projection looks fine as well as the range of the data and it does also overlap with the vector shapefile. I can do the calculations in Qgis using "Zonal Statistics", but I need to do the same calculations on a bunch of raster data, so it is too timeconsuming to run through all the data one at a time with Qgis.

I can open the AET data with GDAL and it looks like this, without any problems:

AET_raster = gdal.Open(filepath)

I thought it was something with the dimensions of the file but this also looks fine. The string of the raster is like so 

<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x119250810> >

And the dimesions are: 

RasterCount = {int} 1, RasterXSize = {int} 1115, RasterYSize = {int} 834

These are the same dimensions a the raster file I have for precipitation.

Does anyone know what can be wrong here?

Idalh
  • 31
  • 3

1 Answers1

0

I had the same issue and fixed it like this:

(i) I opened the .tif with rasterio, accessed the geotransformation and changed the algrebraic signs of "yres" and "ymin":

OLD: 'transform': Affine(0.08333333333334281, 0.0, -180.0,
   0.0, 0.0833333333333286, -90.0)

NEW: 'transform': Affine(0.08333333333334281, 0.0, -180.0,
   0.0, -0.0833333333333286, 90.0)

That would look like in my case:

old_tif = rasterio.open('old.tif')
print(old_tif.profile) # copy & paste the output and change signs 
new_tif_profile = old_tif.profile
new_tif_profile['transform'] = Affine(0.08333333333334281, 0.0, -180.0,
       0.0, -0.0833333333333286, 90.0)

For the sake of simplicity, I copied the old 'transform' and changed the signs manually, you'd do it more correctly by accessing each entry and multiply it with the corresponding sign instead. One option would be with xarray. For knowing which entries to modify, you can have a look at your precipitation data, which seems to work with zonal_stats.

(ii) I then flipped the numpy array over both axis:

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

(iii) Finally, I wrote the file:

with rasterio.open('new.tif', "w", **new_tif_profile) as dest:
dest.write(new_tif_array, indexes=1)

EDIT: The problem is probably due to a "presumption of orthogonality with respect to the Cartesian coordinate axes" within zonal_stats, see Compute zonal statistics for a numpy array in rotated-coordinate space or here: https://github.com/perrygeo/python-rasterstats/issues/98

ConZZito
  • 325
  • 2
  • 5