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.
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:
- create a shapefile containing just the area of interest clipping polygon
ogrinfoto determine the extent of the clipping shapefile
gdal_translateto clip to the shape extents
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.
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, craster.extent craster.fillrasterpoints(xres, yres) craster.getareaofinterest("DATA/area_of_interest.shp") minX, maxX=craster.new_extent -5,craster.new_extent+5 minY, maxY= craster.new_extent -5,craster.new_extent+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...).
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).