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 yOrigin = transform pixelWidth = transform pixelHeight = transform # 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.