Today I had to clip a raster image by 102,000 polygon features. With out doubt this task demanded automation. I decided to try my hand at python coding using the excellent GDAL library. I had this great idea that i could do the whole thing nicely in GDAL, but i was foiled by time (or not enough of it) and I did half of the process in a pure GDAL script then I dropped to the cmd line to execute gdal_translate with a changing projwin per iteration. In case you are unaware with the gdal library comes a series of astonishingly useful utilities. I used to use the FWTools port but now I use the OSGeo4w one instead. There are a number of raster processing wonder-tools in here, and gdal_translate is one of them! Back to my script, I realize now that in actual fact i am using almost exclusively OGR pieces of GDAL, confused yet?. Its simple: GDAL = rasters OGR = Vectors I start by importing:

import osgeo.gdal as gdal
import osgeo.ogr as ogr
import os

Then I open up my database polygon layer and print off the extents and the number of features so I’m happy the layer is good:

driver = ogr.GetDriverByName('MYSQL')
vec = "MYSQL:database,user=root,tables=my_tables"
vector_ds = driver.Open(vec, 0)
layer = vector_ds.GetLayer()		
numFeatures = layer.GetFeatureCount()
print 'Feature count: ' + str(numFeatures)
extent = layer.GetExtent()
print 'Extent:', extent
print 'UL:', extent[0], extent[3]
print 'LR:', extent[1], extent[2]

feature = layer.GetNextFeature()
count = 0
while feature:
	feature_id = feature.GetFieldAsString('objectid')
	geometry = feature.GetGeometryRef()

	coords = str(geometry)
	coords1 = coords[10:-2].split(",")
	co_dic = {}
	count2 = 1
	
	ulx = 0
	uly = 0
	lrx = 0
	lry = 0
	
	for c in coords1:
		d = c.split(" ")
		co_dic = {count2 : {'X': float(d[0]), "Y":float(d[1])}}
		
		#print co_dic
		co_dic
		
		if count2 == 1:
			ulx = co_dic[count2]["X"]
			uly = co_dic[count2]["Y"]
			lrx = co_dic[count2]["X"]
			lry = co_dic[count2]["Y"]
			
		if co_dic[count2]["X"] < ulx:
			ulx = co_dic[count2]["X"]
		
		if co_dic[count2]["X"] > lrx:
			lrx = co_dic[count2]["X"]
		
		if co_dic[count2]["Y"] > lry:
			uly = co_dic[count2]["Y"]
			
		if co_dic[count2]["Y"] < lry:
			lry = co_dic[count2]["Y"]
			
			
		count2 += 1
	coords2 = coords1[0].split(" ")[0]
	
	print "*** Processing Parcel " + feature_id  + " ***"
	print "*** ulx= " + str(ulx) + " :: uly= " + str(uly) + " :: lrx= " + str(lrx) + " :: lry= "  + str(lry)
	cmd = "gdal_translate -projwin " + str(ulx) + " " + str(uly) + " "  + str(lrx) + " "  + str(lry) + " C:/image.tif " + feature_id + ".tif"
	
	os.popen(cmd)

	feature = layer.GetNextFeature()
	count = count + 1