Weighted zonal statistics in Python

Question

I'm asking something similar than Weighted zonal statistics in Python and/or R. I know very well the R-side of this question, but I can't figure out the answer for Python.

So far I can compute non-weighted zonal statistics:

import os
import rioxarray as rioxr
from rioxarray.merge import merge_arrays
from glob import glob
import geopandas

polys = geopandas.read_file(path_to_polygons)

files = glob(os.path.join('path/to/MOD13Q1/2002/033','*.hdf')) # example path
nbr_single = []
for file in files:
    with rioxr.open_rasterio(file) as src:
        NIR = getattr(src,'250m 16 days NIR reflectance')
        MIR = getattr(src,'250m 16 days MIR reflectance')
        nbr_single.append((NIR-MIR)/(NIR+MIR))

nbr_mosaic = merge_arrays(nbr_single) 

nbr_zonal = zonal_stats(
    vectors=polys,
    raster=nbr_mosaic[0,:,:].values, # since it's a 3D array
    affine=nbr_mosaic.rio.transform(),
    stats=["sum"],
)

# for testing in R
nbr_mosaic.rio.to_raster(os.path.join(some_path,'mosaic.tif'), compress='LZW')

And I tried a weighted and non-weighted approach in R:

library(terra)
library(sf)
library(exactextractr)

r <- rast('some_path/to/mosaic.tif')
sf <- read_sf(path_to_polygons)
v_ <- vect(path_to_polygons)

r1_values <- exact_extract(r,sf, fun = 'count')
r2_values <- extract(r,v_,length)

Where r1_values is a weigthed approach and r2_values is a non-weighted (although terra::extract can compute weighted statistic as well)

R's weighted count vs Python's non-weighted count:

R's non-weighted count vs Python's non-weighted count (minimal difference):

How can I compute a weight raster or directly weigthed zonal statistics in python?

Answer 1

Now exactextract is available for Python:

pip install exactextract

Basic operation example:

import geopandas as gpd
import rioxarray as rxr
import exactextract

gdf = gpd.read_file('/path/to/vector.shp')
r = rxr.open_rasterio('/path/to/raster.tif')
dict_result = exactextract.exact_extract(r, gdf, 'mean')