A simpler solution, although less elegant is simply use rasterio, gdal and numpy. Let's say that in our shapefile we have a column called 'column' containing the features names (our polygons), let's say that our polygons are 'A', 'B' and 'C'.
import gdal, os, rasterio
import numpy as np
im = 'image.tif'
shp = 'shapefile.shp'
#open the image using rasterio and reshape it to a np.array
a = rasterio.open(im).read()
a = a.reshape((a.shape,a.shape))
# gdal is weird when operating with folders sometimes,
# so we just move to the folder where the image is
os.chdir(folder) #folder with image
#now we run gdalwarp with the appropriate parameters in a loop
# over our polygons. This is the part that could be improved
# by using fiona or geopandas, something to extract the features
# names and dump them in a list.
# The trick with the features names is to use single AND double quotes, like "'this'", because gdal wants to get something like 'this' as a parameter.
polnames = ["'A'","'B'","'C'"]
for p in polnames:
os.system('gdalwarp -cutline %s -cwhere "%s = %s" -crop_to_cutline %s %s_%s' %(shp,c,p,im,p,im))
Now you have the image cut in pieces according to your features in the shapefile. Then, we just open the images and get the data.
# open one image using rasterio
b = rasterio.open('Aim.tif').read().reshape((b.shape,b.shape))
# if you want, you can flatten and have a 1d array
b = b.flatten()
Done! I said it is not elegant, there is a lot of room to improve, but it works.