I have downloaded multiple products from EO-1 satellite and I have created digital rasters out of them. Each of them contains 10 or 242 bands (those are hyperspectral). I have also created a shapefile for each of them with multiple points. Now I would like to create a table/array with Python containing values of a pixel from each band with coordinates from a shapefile. I would like to have all bands in separate columns and all points in separate rows. I have managed to create such data frame in R, and I tried to do so with GDAL and numpy but I am far less familiar with Python.