0

I am trying to apply the rasterio.features.rasterize function in order to obtain a zonal statistics from a vector file containing multiple vectors.

Despite my efforts, I have been unsuccessful in this task. All I can do is rasterize one single shapefile-feature-file per time.

Would anyone know how to properly rasterize multiple features using rasterio?

My ultimate goal would be an algorithm so that:

For each iteration, each vector feature would produce a single matriz mask in memory. This memory-matriz would then be applied over my raster data with same dimensions as that one. Zonal statistics as (mean, maximum, minimum, std, var) would then be inserted as new properties of each rasterized feature of my original vector (shapefile) file.

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
Philipe Riskalla Leal
  • 179
  • 7
  • 3
    Can you show how you do it for one feature per time? Once you are there, it should not be to hard to automatically do it for all features of your shapefile – joris Feb 19 '18 at 14:55
  • with your current approach you probably want to loop through the features in your shapefile. – RutgerH Feb 19 '18 at 15:27

1 Answers1

1

Based on Extract raster values within shapefile with pygeoprocessing or gdal, we can write the following function to extract values from the raster dataset based on a geometry:

from shapely.geometry import mapping
from rasterio.mask import mask

def extract_zonal_statistics(src, geom, stat, band):
    geom_map = mapping(geom)
    out_image, out_transform = mask(src, [geom_map], crop=True, filled=False)
    return stat(out_image[band - 1])

Considering gdf being our GeoDataFrame, we can then add a column with zonal statistics:

import rasterio

with rasterio.open("/path/to/raster.tif") as src:
    gdf['zonal_stat_mean'] = gdf.geometry.apply(
        lambda geom: extract_zonal_statistics(src, geom, stat=np.mean, band=1))

Alternatively, you could use the rasterstats library which also implements this functionality, and combine the results with the original frame:

from rasterstats import zonal_stats

with rasterio.open("/path/to/raster.tif") as src:
    affine = src.transform
    array = src.read(1)
    df_zonal_stats = pd.DataFrame(zonal_stats(gdf, array, affine=affine))

gdf2 = pd.concat([gdf, df_zonal_stats], axis=1)
joris
  • 3,903
  • 22
  • 28
  • Dear Joris,thank you for you reply. I have tried this code on my computer, but an error appeared: """ File "C:\Anaconda3\lib\site-packages\rasterstats\io.py", line 134, in rowcol r = int(op((y - affine.f) / affine.e))

    ZeroDivisionError: float division by zero"""

    Despite my efforts in verifing where this division occured in the main code, I couldn't solve it. Would you perhaps know what this problem is, and how to solve it? Sincerely, Philipe Leal

     – Philipe Riskalla Leal Feb 20 '18 at 20:30
  • Dear Joris, I fortgot to mention that the first script did work just fine, once I changed the extract_zonal_statistics function to:
    def extract_zonal_statistics(src, geom, stat, band): geom_map = mapping(geom) out_image, out_transform = mask(src, [geom_map], crop=True, nodata=np.nan, all_touched=True) return stat(out_image[band - 1])     – Philipe Riskalla Leal Feb 20 '18 at 20:48