Determining if shapefile and raster overlap in Python using OGR/GDAL?
The following script determines the bounding box of a raster and creates based on the bounding box a geometry.
import ogr, gdal
raster = gdal.Open('sample.tif')
vector = ogr.Open('sample.shp')
# Get raster geometry
transform = raster.GetGeoTransform()
pixelWidth = transform[1]
pixelHeight = transform[5]
cols = raster.RasterXSize
rows = raster.RasterYSize
xLeft = transform[0]
yTop = transform[3]
xRight = xLeft+cols*pixelWidth
yBottom = yTop-rows*pixelHeight
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(xLeft, yTop)
ring.AddPoint(xLeft, yBottom)
ring.AddPoint(xRight, yTop)
ring.AddPoint(xRight, yBottom)
ring.AddPoint(xLeft, yTop)
rasterGeometry = ogr.Geometry(ogr.wkbPolygon)
rasterGeometry.AddGeometry(ring)
Next, the geometry of the vector polygon is determined. This answers your first question.
# Get vector geometry
layer = vector.GetLayer()
feature = layer.GetFeature(0)
vectorGeometry = feature.GetGeometryRef()
Last, the geometry of the vector and raster are tested for intersection (returns True
or False
). This answers your second question.
print rasterGeometry.Intersect(vectorGeometry)
I find @ustroetz solution very helpful but it needed to be corrected in two places. Firstly, pixelHeight = transform[5] is already negative value, so equation should be
yBottom = yTop+rows*pixelHeight
Secondly, the order of the points in the ring must be counter clockwise. I was having problems with that. Correct order is:
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(xLeft, yTop)
ring.AddPoint(xLeft, yBottom)
ring.AddPoint(xRight, yBottom)
ring.AddPoint(xRight, yTop)
ring.AddPoint(xLeft, yTop)