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]

Then I iterate through the features and do some nasty hackish string processing to get my coordinates working ok:

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
		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"

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

Did you catch that third last line… yeah, there’s the hack, sorry folks. The lot of the consultant is: “you do your best when you have time, but when it runs out you just make it work” Its 3 hours later and I’m on feature 89561. I love python!

Sharing options
mountain during golden hour

Observations • Will Cadell

Sparkgeo & N51

Sparkgeo is very excited to be co-hosting N51 in Banff this year. Come and let the scenery take your breath away, while N51 stimulates your geospatial ideas.

selective focus photography of licensed plate with open text hanged

Observations • Dan Ormsby

FOSS4G UK 2022 – Sparkgeo represent at key open-source conference

The UK “free and open-source software for geo” (FOSS4G) community came together on November 17th to celebrate PostGIS Day.

closeup photography of clear glass window closed

Observations • Darren Wiens

How to use SVG Filters on Web Maps

Here at Sparkgeo, while we often prefer to provide answers over visualizations, we still make a lot of web maps. And when we do, we take pride in…

Need a geospatial partner?

Our team complements organizations like yours—by providing on-tap access to geospatial, analytics, and mapping expertise.

Let’s talk

Join our team?

We’re always looking for skilled technologists to help us build the future of geospatial. Got a minute to find out more about us?

Working Here

Sharing options