0

I'm trying to use the GeometryReference.Intersects() method to see if a feature within a grid intersects with another polygon. I have tried it in the QGIS python console (windows) and directly on linux, and it crashes in both cases.

EDIT: I have tried the same code with the data reprojected to UTM, which is working. So possibly it is a problem with using the Intersect method the the OSGB projection (EPSG: 27700)

If possible I would like to know what is going wrong, or alternatively if there is another method of checking if two features intersect within python.

Here is the code I am trying:

from osgeo import ogr

driver = ogr.GetDriverByName("ESRI Shapefile")

grid_ds = driver.Open("test_grid.shp")

poly_ds = driver.Open("test_poly.shp")

grid_layer = grid_ds.GetLayer()
poly_layer = poly_ds.GetLayer()

poly_geom = poly_layer[0].GetGeometryRef()
grid_geom0 = grid_layer[0].GetGeometryref()

grid_geom0.Intersects(poly_geom)

Here is a link to some example datasets

1

You need to hold a reference to the feature you are grabbing. It's described more here:

http://trac.osgeo.org/gdal/wiki/PythonGotchas

from osgeo import ogr

driver = ogr.GetDriverByName("ESRI Shapefile")

grid_ds = driver.Open("test_grid.shp")

poly_ds = driver.Open("test_poly.shp")

grid_layer = grid_ds.GetLayer()
poly_layer = poly_ds.GetLayer()

poly_feat = poly_layer[0]
grid_feat = grid_layer[0]

poly_geom = poly_feat.GetGeometryRef()
grid_geom0 = grid_feat.GetGeometryRef()

print(grid_geom0.Intersects(poly_geom))

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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