Grid Shape & Multiple TIF

[manuvare]

First at all, many thanks on sharing all this work.

I´m new on Python and maybe i want to run before walk.

I´m running the script orienting to climate data, where I have a Shape with around 500 equal squares and 30 years daily data of rainfall (tif files for each day around 11k .tif).

I need as an output a table with the rainfall value for each day for each grid.

How to deal with this?

Thanks

[perrygeo]

I would loop over the climate tifs and run zonal stats on each:

rasters = [
    'rainfall_001.tif',
    'rainfall_002.tif'
]

results = {}

for raster in rasters:
    results[raster] = zonal_stats(shapefile, raster, stats=['mean'])

You would be left with a results dictionary like

{
    'rainfall_001.tif': [     # a list of stats for each raster
        {'mean': 123.44},  # one stats record per feature in shapefile
        {'mean': 234.55},
        ...
    ]
...
}

From here, how you process this data is entirely up to your application.

[manuvare]

I think i have something wrong with packages:

FutureWarning: The value of this property will change in version 1.0. Please see rasterio/rasterio#86 for details.
with Raster(raster, affine, nodata, band) as rast:

FutureWarning: GDAL-style transforms are deprecated and will not be supported in Rasterio 1.0.
self.affine = guard_transform(self.src.transform)

UserWarning: Setting nodata to -999; specify nodata explicitly
warnings.warn("Setting nodata to -999; specify nodata explicitly")

Im getting results. However, im getting different values that the ones i could get using qgis on some (not all) grids .

[perrygeo]

The transforms warnings should not cause any problems, it's a known issue for future releases #64

The user warning re: nodata indicates that you may need to set nodata explicitly.

im getting different values that the ones i could get using qgis

Can't help you unless you provide more information. First thing to check is making sure you're set the nodata correctly. Please open a new issue if you can find and reproduce a bug.

[manuvare]

It was a double check but different are not sensible at this stage.

Nodat set now is ok.

I will check how to include on rasters tuple all the files in my directory (up to 10k) and how to export results to a csv.

Thanks!

[manuvare]

Done.

IS there any way to get the code more efficient?

import os, glob
import rasterstats, rasterio
from rasterstats import zonal_stats

polys = r'C:\Users\mvarela3\Documents\CHIRPS\Python\GRIDBSAS.shp'

rasterfolder = r'C:\Users\mvarela3\Documents\CHIRPS\img\Global'
os.chdir(rasterfolder)

results = {}

for lyr in glob.glob("*.tif") :
results[lyr] = zonal_stats(polys, lyr, stats=['mean'], nodata=-999,all_touched=True, geojson_out=True)

with open('my_file.csv', 'w') as f:
[f.write('{0},{1}\n'.format(key, value)) for key, value in results.items()]