zonal_stats does not include expected raster cells

[Tinkaa]

I created a toy xr.DataArray and gpd.GeoDataFrame. When applying zonal_stats to this, the results are not as expected. It identifies 2 cells within one polygon and 1 within the other. I would expect the result to have 4 cells in the one polygon, and 0 in the other. I have tried to dig into the code what causes this but haven't been able to understand it. I am guessing it has something to do with the affine or how the polygons are defined. When I apply rio.clip, it does work as expected. Do you know what causes this unexpected behaviour of zonal_stats?

All code is shown below

import xarray as xr
from shapely.geometry import Polygon
import geopandas as gpd
import matplotlib.pyplot as plt
from rasterstats import zonal_stats
import rioxarray

da = xr.DataArray(
[[1,2,3],
       [4,5,6]],
dims=("x", "y"),
coords={"x": [0.5,1.5], "y": [2.5,1.5,0.5]},
)

d = {'geometry': [Polygon([(0, 0), (0,2), (2, 2), (2, 0)]),Polygon([(2,0),(2,2),(3,2),(3,0)])]}
gdf = gpd.GeoDataFrame(d)

zonal_stats(
    vectors=gdf,
    raster=da.values,
    affine=da.rio.transform(),
    stats=["mean","count"],
    all_touched=True,
)

Outputs: [{'mean': 4.5, 'count': 2}, {'mean': 6.0, 'count': 1}]

#plotting to illustrate the data
fig, ax = plt.subplots()
da.plot(ax=ax,x="x",y="y")
gdf.boundary.plot(ax=ax, color='r', lw=1)
ax.set_xlim(0,4)
ax.set_ylim(0,4)