9

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?

closed as too broad by PolyGeo Jul 15 '17 at 1:40

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer. Avoid asking multiple distinct questions at once. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.

  • 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
15

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)
  • 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
3

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)

Not the answer you're looking for? Browse other questions tagged or ask your own question.