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, systry: 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 layerdriver = 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 driversgdal.AllRegister()# open the imageimg = gdal.Open('b5.TIF', GA_ReadOnly)if img is None: print 'Could not open landsat.TIF' sys.exit(1)# get image sizerows = img.RasterYSizecols = img.RasterXSizebands = img.RasterCount# get georeference infotransform = img.GetGeoTransform()xOrigin = transform[0]yOrigin = transform[3]pixelWidth = transform[1]pixelHeight = transform[5]# loop through the features in the shapefilefeat = 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 shapefileshp.Destroy()
أكثر...
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, systry: 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 layerdriver = 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 driversgdal.AllRegister()# open the imageimg = gdal.Open('b5.TIF', GA_ReadOnly)if img is None: print 'Could not open landsat.TIF' sys.exit(1)# get image sizerows = img.RasterYSizecols = img.RasterXSizebands = img.RasterCount# get georeference infotransform = img.GetGeoTransform()xOrigin = transform[0]yOrigin = transform[3]pixelWidth = transform[1]pixelHeight = transform[5]# loop through the features in the shapefilefeat = 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 shapefileshp.Destroy()
أكثر...