Clipping polygon precisely to raster (extent) using GDAL?

Question

I need to clip a vector layer to the extent of a raster. It has to be very precise that's why I would not want to draw a polygon myself (and clip the two polygons). There seem to be some ways to clip a raster to a polygon but I have difficulties to find a solution to clip a polygon to a raster. I already tried to create a polygon out of my raster several times in ArcGIS and R, so I could clip the result. But every time it crashes because the raster is to big. There must be a simple way to solve this, and somehow I can just not see it at the moment.

How can I do this using GDAL?

Answer 1

You can use GDAL's ogr2ogr utility (http://www.gdal.org/ogr2ogr.html) to clip a polygon by an extent. If you have GDAL for Python you can use something similar to this script to perform the task:

from osgeo import gdal
import subprocess

# define file paths
inRasterPath = 'C:/path/to/file.tif'
inVectorPath = 'C:/path/to/file.shp'
outVectorPath = 'C:/path/to/file.shp' 

# get the extent of the raster
src = gdal.Open(inRasterPath)
ulx, xres, xskew, uly, yskew, yres  = src.GetGeoTransform()
sizeX = src.RasterXSize * xres
sizeY = src.RasterYSize * yres
lrx = ulx + sizeX
lry = uly + sizeY
src = None

# format the extent coords
extent = '{0} {1} {2} {3}'.format(ulx, lry, lrx, uly);

# make clip command with ogr2ogr - default to shapefile format
cmd = 'ogr2ogr ' + outVectorPath + ' ' + inVectorPath + ' -clipsrc ' + extent

# call the command
subprocess.call(cmd, shell=True)

If you don't have GDAL for Python, you can use the OSGeo4W shell that comes with QGIS and run the gdalinfo command on your raster to get the extent and then the ogr2ogr command to clip the vector to the extent noted from the gdalinfo result with the -clipsrc flag.

Make sure your polygon and raster file are in the same projection.