12

I am building a script in python using OGR/GDAL.

I have a set of shapefiles and a set of GeoTiff raster files.

I would like to have my script ignore shapefiles if they do not intersect with the raster area.

The shapefile is not a rectangle, so I cannot simply compare the xmin/xmax,ymin/ymax values returned by layer.GetExtent(). I need the actual polygon representing it's overall shape, and then some way of determining if that polygon intersects with the raster square.

I was thinking I could somehow merge all the polygons in the shapefile into one feature, and then read the geometry on that feature, and then compare that information to the raster extent. However, I am unsure of specifically how to execute this.

  1. How to extract border polygon information from shapefile?
  2. How to determine if that polygon intersects a given square area?
PolyGeo
  • 65,136
  • 29
  • 109
  • 338
JFerg
  • 347
  • 1
  • 2
  • 9
  • I'm not familiar with osgeo, but the arcpy equivalent would (might) involve: read raster extent, create polygon covering extent in memory, cycle through shapefiles, clipping each to extent rectangle, test if anything results. – phloem Dec 12 '14 at 21:59

2 Answers2

20

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)
ustroetz
  • 7,994
  • 10
  • 72
  • 118
  • 2
    Thanks, this was exactly what I was looking for. This answer clearly shows how to create, extract, and run functions between geometry objects, which is exactly what I was looking for. – JFerg Dec 15 '14 at 15:01
  • This solution test whether the FID = 0 polygon intersects with the raster. How do you get the geometry of the shapefile total when no polygon represents this? – JFerg Jan 29 '15 at 19:51
  • 1
    EDIT: The increase in computation time is inconsequential so I check if each polygon in the shapefile intersects now. – JFerg Jan 29 '15 at 21:59
  • Thanks for the helpful solution, however pixelHeight seems to be negative when the zone is in the Southern hemisphere (we should check the behavior of raster.GetGeoTransform()). yBottom = yTop-abs(rows*pixelHeight) fixes this. – Franck Theeten Apr 02 '22 at 14:35
6

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)
Pandza
  • 163
  • 1
  • 7