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

  from osgeo import ogr, gdal
  from osgeo.gdalconst import *
except ImportError:
  import ogr, gdal
  from gdalconst import *

# 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'
shpLayer = shp.GetLayer()

# register all of the GDAL drivers

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

# 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 = shpLayer.GetNextFeature()

# close the shapefile

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.

closed as off-topic by PolyGeo Nov 7 '15 at 22:33

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "Questions seeking help to debug/write/improve code must include the desired behavior, a specific problem or error and the shortest code necessary to reproduce it in the question itself. Providing a clear problem statement and evidence of a code attempt will help others to help you. See: How to create a Minimal, Complete, and Verifiable example." – PolyGeo
If this question can be reworded to fit the rules in the help center, please edit the question.


The title ('Extract raster values ...') is not reflecting your code instructions because you are mixing vector and raster layers. Assuming that the answer is only for raster layers, I think that the best way to read in them (for any purpose with Python/GDAL) is by using a scanline and the unpack struct function. The code is more compact, the control is more effective and the execution time is faster than the one with 'ReadAsArray'. The scanline/struct method depends on fmttypes and their values can be supplied in a dictionary. An example of code (searching by rows) can be observed below:

from osgeo import gdal
import os
import struct

fmttypes = {'Byte':'B', 'UInt16':'H', 'Int16':'h', 'UInt32':'I', 'Int32':'i', 'Float32':'f', 'Float64':'d'}

path= "/home/zeito/Desktop/CURSO_QGIS/aleatorio.tif"

(raiz, filename) = os.path.split(path)

dataset = gdal.Open(path)

band = dataset.GetRasterBand(1)

print "rows = %d columns = %d" % (band.YSize, band.XSize)

BandType = gdal.GetDataTypeName(band.DataType)

print "Data type = ", BandType

print "Processing %s" % filename
print "whose path is %s" % raiz

for y in range(band.YSize):

    scanline = band.ReadRaster(0, y, band.XSize, 1, band.XSize, 1, band.DataType)
    values = struct.unpack(fmttypes[BandType] * band.XSize, scanline)
    print values

dataset = None

When this code was run using a litle raster of 29 rows by 29 columns, for a better visualization, I got the next result (in Windows you can use the Python IDLE):

enter image description here

I only write the values at the screen but they can be used for being modified and create another raster; for example.

Not the answer you're looking for? Browse other questions tagged or ask your own question.