Getting negative values after a zonal_stats in Python

Question

I have a shapefile with overlapping polygons and raster data about suitability, both of them have the same CRS.

I would like to summarize some statistics on my polygons, using the "all touched" rasterization strategy. However, for some polygons I obtain a negative value even although the raster values should be greater or equal than 0. When I repeat the zonal statistics without "all touched" option, those negative values either become 0 or missing (.).

I would like to understand why the zonal stats assign negative values to some polygons choosing the "all touched" rasterization strategy, but without it 0 or missing values.

This is an example of the code I am using:

import numpy as np
import pandas as pd
import geopandas as gpd
from pyproj import CRS
import rasterio
from rasterstats import zonal_stats
from rasterio.mask import mask

polygons = gpd.read_file(data_directory + 'GridCells/ccode_GID.shp')
raster = rasterio.open(tif_file)

raster.profile
{'driver': 'GTiff', 'dtype': 'int16', 'nodata': -9.0, 'width': 4320, 'height': 2160, 'count': 1, 'crs': CRS.from_epsg(4326), 'transform': Affine(0.08333333333333333, 0.0, -180.0,
       0.0, -0.08333333333333333, 90.0), 'blockxsize': 256, 'blockysize': 256, 'tiled': True, 'compress': 'lzw', 'interleave': 'band'}

polygons.crs

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich


raster_array, raster_transform = mask(raster, shapes=polygons.geometry, crop=True, nodata=np.nan)

zs1 = zonal_stats(polygons, raster_array[0], affine=raster_transform, stats=['mean'], all_touched=True)
polygons['value'] = [x['mean'] for x in zs1]