GDAL Clipping (hack)

Question

I am trying to use the python script from http://www.sparkgeo.com/labs/2010/07/07/gdal-clipping-hack/ to clip a raster based on many shapefiles. The script is using a mysql database an i would like to use it on a directory of shapefiles. I am not a programmer so I'am not sure how to make this work.

Answer 1

import os
from osgeo import gdal
from osgeo import ogr    
os.chdir('C:/yourpath')
files = glob.glob('*.shp')
for item in files:
    ''' This will loop over all shapefiles in the directory specified'''
    driver = ogr.GetDriverByName('ESRI Shapefile')
    vector_ds = driver.Open(item, 0)
    layer = vector_ds.GetLayer()
    #...
    #At the end of the loop do the following to clean up
    vector_ds.Destroy()

Then you can pick up with the rest of the script in that post, putting it into the ...

Not knowing 100% what you want to do keep in mind the QGIS GDALTools - Clipper will show you how to build a gdalwarp command to do something similar.

Answer 2

Wow! That's an incredibly inefficient way of clipping a raster. OGR/GDAL can be used straight from the command line without recourse to Python in this case.

First, you need to merge your shapefiles into one big one. Assuming you're using Windows, the following DOS batch file is taken from here (a UNIX/Linux script will be syntactically different, but will follow the same form):

mkdir merged
for %f in (*.shp) do (
if not exist merged\merged.shp (
ogr2ogr -f "esri shapefile" merged\merged.shp %f) else (
ogr2ogr -f "esri shapefile" -update -append merged\merged.shp %f -nln Merged )
)

This generates a shapefile called merged.shp with a single layer called Merged.

Then you can clip your raster based on this shapefile quite simply using gdalwarp like this:

gdalwarp -cutline merged\merged.shp -cl Merged <source image> <dest image>