Friday 17 June 2016

shapefile - Clipping polygon precisely to raster (extent) using GDAL?


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




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.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...