python - Extract raster values gdal looping rasters

I have this code, most of them getting from this forum. I need to loop over a folder and sub-folders, or just to run a list (putting previously all the raster inside, I don't mind one way or the other).

But I don't find out how to do the loop, I have done it with arcpy in just 30 minutes, but I am fighting with gdal for a couple of days and nothing...

import os, sys

try:
  from osgeo import ogr, gdal
  from osgeo.gdalconst import *
  os.chdir('/home/digd/Desktop/puntos')
except ImportError:

  import ogr, gdal
  from gdalconst import *
  os.chdir('/home/digd/Desktop/puntos')

# open the shapefile and get the layer
driver = ogr.GetDriverByName('ESRI Shapefile')
shp = driver.Open('/home/digd/Desktop/puntos/centroids2.shp')
if shp is None:
  print 'Could not open puntos.shp'
  sys.exit(1)

shpLayer = shp.GetLayer()

# register all of the GDAL drivers
gdal.AllRegister()



# open the image
img = gdal.Open('b5.TIF', GA_ReadOnly)
if img is None:

  print 'Could not open landsat.TIF'
  sys.exit(1)

# get image size
rows = img.RasterYSize
cols = img.RasterXSize
bands = img.RasterCount

# get georeference info
transform = img.GetGeoTransform()

xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = transform[5]

# loop through the features in the shapefile
feat = shpLayer.GetNextFeature()
while feat:
  geom = feat.GetGeometryRef()
  x = geom.GetX()

  y = geom.GetY()

  # compute pixel offset
  xOffset = int((x - xOrigin) / pixelWidth)
  yOffset = int((y - yOrigin) / pixelHeight)

  # create a string to print out
  s = feat.GetFieldAsString('Name') + ' '

  # loop through the bands

  for j in range(bands):
    band = img.GetRasterBand(j+1) # 1-based index

    # read data and add the value to the string
    data = band.ReadAsArray(xOffset, yOffset, 1, 1)
    value = data[0,0]
    #print value
    s = s + str(value) + ' '

  # print out the data string

  print s

  # get the next feature
  feat.Destroy()
  feat = shpLayer.GetNextFeature()

# close the shapefile
shp.Destroy()

---

Uhmm no, I was trying to do a better explanation of my question.

Anyway, I eventually could do it using numpy and GDAL. I needed to get the pixels values within some polygons, so I convert the polygons to rasters and then I could easily perform different masks and calcs.

I find this process really easy with arcpy, but no idea of how to do this with other tools (well, it's easier and faster and better (at least for me) with the rsgislib, but I have to work in a window enviroment and it's not easy to install this library in windows (at least for me).

Again, I don't know of this is an answer or whatever it's (and truly I don't care), you asked me something and I am just trying to give you the best answer that I can.