Clipping raster with vector layer using GDAL

Question

I have installed GDAL using the Osgeo installer. How can I clip a raster layer with a vector layer programmatically? Is there any GDAL API available which can help me with this? I am using Python.

Answer 1

I'm not sure about the gdal api, there is void* GDALWarpOptions::hCutline in the Warp Options referenced from the Warp API tutorial, but no explicit examples. Are you sure you need a programmatic answer? The command line utilities can do it out of the box:

1. create a shapefile containing just the area of interest clipping polygon
2. use ogrinfo to determine the extent of the clipping shapefile
3. use gdal_translate to clip to the shape extents
4. use gdalwarp with -cutline parameter

Steps 2 & 3 are for optimization, you could get by with just gdalwarp -cutline ....

See Clipping rasters with GDAL using polygons from Linfinity for linux based solution all wrapped up in one script. Another cutline example can be seen in Michael Corey's tutorial creating hillshades for Mapnik.

Answer 2

Joel Lawhead of GeospatialPython has complete python example in Clip raster using shapefile, a well written tutorial. You'll need to install the Python Image Library (PIL) which isn't included in Osgeo4W (for which you might need to add o4w python to windows registry to get the install program to work).

Answer 3

It seems that this subject is always coming back. I myself didn't know that GDAL >1.8 is so advanced it already gives you fair command line handling to do that task.

The comment from Mike Toews is pretty useful but you could simply do for example:

gdalwarp -of GTiff -cutline DATA/area_of_interest.shp -cl area_of_interest  -crop_to_cutline DATA/PCE_in_gw.asc  data_masked7.tiff 

You could wrap this command inside a python script with the excellent subprocess module.

One thing which was really problematic for me is that I needed to supply a minimal solution to that problem, meaning as simple as possible and does not require to many external dependecies. The use of Python Imaging Library as in the tutorial by Joel Lawhead is neat, but I came up with the following solution: using Numpy masked arrays.
I don't know if it's better, but that was what I knew back than (3 years ago...).
Originally I created a valid data area inside the original raster (e.g. the extent of the output raster where the same), but I liked the idea of making the raster also smaller (e.g. -crop_to_cutline), so I adopted world2Pixel from Joel Lawhead. Here is my own solution:

def RasterClipper():
    craster = MaskRaster()
    contraster2 = 'PCE_in_gw.aux'
    craster.reader("DATA/"+contraster2.replace('aux','asc'))
    xres, yres = craster.extent[1], craster.extent[1]
    craster.fillrasterpoints(xres, yres)
    craster.getareaofinterest("DATA/area_of_interest.shp")
    minX, maxX=craster.new_extent [0]-5,craster.new_extent[1]+5
    minY, maxY= craster.new_extent [2]-5,craster.new_extent[3]+5
    ulX, ulY=world2Pixel(craster.extent, minX, maxY)
    lrX, lrY=world2Pixel(craster.extent, maxX, minY)
    craster.getmask(craster.corners)
    craster.mask=np.logical_not(craster.mask)
    craster.mask.resize(craster.Yrange.size,craster.Xrange.size)
    # choose all data points inside the square boundaries of the AOI,
    # replace all other points with NULL
    craster.cdata= np.choose(np.flipud(craster.mask), (craster.data, -9999))
    # resise the data set to be the size of the squared polygon
    craster.ccdata=craster.cdata[ulY:lrY, ulX:lrX]
    craster.writer("ccdata2m.asc",craster.ccdata, (minX+xres*.5, maxY+yres*.5), 10,10,Flip=False)
    # in second step we rechoose all the data points which are inside the
    # bounding vertices of AOI
    # need to re-define our raster points
    craster.xllcorner, craster.yllcorner = minX, minY
    craster.xurcorner, craster.yurcorner = maxX, maxY
    craster.fillrasterpoints(10,10)
    craster.getmask(craster.boundingvertices) # just a wrapper around matplotlib.nxutils.points_in_poly
    craster.data=craster.ccdata
    craster.clip2(new_extent_polygon=craster.boundingvertices)
    craster.data = np.ma.MaskedArray(craster.data, mask=craster.mask)
    craster.data = np.ma.filled(craster.data, fill_value=-9999)
    # write the raster to disk
    craster.writer("ccdata2m_clipped.asc",craster.data, (minX+xres*.5, maxY+yres*.5), 10,10,Flip=False)

for a full description of the class MaskRaster and it's methods, see my project's github.

Using this code you will still need to use GDAL. However, the plan is to use in the future pure Python where I can, because the intended audience of my software has difficulties with too many dependencies (I use Debian to develop the software, and the clients use Windows 7...).