My main goal is to get the percentage of vegetation cover of gardens. The figure below shows the gardens (orange borders) with vegetation and non-vegetation.

Snapshot of used data:

I want for each garden the percentage of vegetation cover (number of green pixels /total number pixels. This computation needs to be done in Python. In QGIS I found a tool called 'Zonal histogram' which can be used for this.

Do you know a Python way to achieve the same or maybe how I could create a QGIS script which can be called by Python and runs this QGIS tool with output from Python?


rasterstats zonal_stats with categorical=True.

For example with pandas so you can manipulate the data as you like:

import geopandas as gpd
import pandas as pd
from rasterstats import zonal_stats

lots = gpd.read_file(r'/home/bera/Desktop/tempgis/lots.shp')
raster = r'/home/bera/Desktop/tempgis/S2_maxlik.tif'

zs = zonal_stats(vectors=lots['geometry'], raster=raster, categorical=True)
stats = pd.DataFrame(zs).fillna(0) #One column per raster category, and pixel count as value
stats_percentage = stats.apply(lambda x: round(x/x.sum()*100,1), axis=1) #From pixel count to percentage

#      2     4     1    3     0
#0  99.6   0.4   0.0  0.0   0.0
#1  71.1  13.3   6.2  9.3   0.0
#2  37.8  34.2  28.0  0.0   0.0
#3  73.3   0.0   0.0  0.0  26.7

#Rename columns
stats_percentage = stats_percentage.rename(columns={0:'class_0',
                                                      4:'class_4'}, errors="raise")

#Merge(=join) the results back to your input
results = pd.merge(left=lots, right=stats_percentage, how='left', left_index=True, right_index=True)
results = results.reindex(sorted(results.columns), axis=1) #Sort columns by name

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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