python - Extract raster values within shapefile with pygeoprocessing or gdal

I would like to know how to get all the raster values within a polygon using gdal or pygeoprocessing, without reading the entire grid as an array.

pygeoprocessing and gdal can do zonal statistics but only the min, max, mean, stdev or count are available from such a function. Since zonal statistics need to access the values, would it be easy to extract values the same way ?

I found a very similar question here : (Getting pixel value of GDAL raster under OGR point without NumPy?) but only for a particular "point".