GDAL: Get Feature Names of intersecting Shp-Files

Question

I am looking for all Featurenames in a shapefile (shp2) which intersect with a given Featurename ('MyFeatureName') in another shapefile (shp1).

The following code works (python), however I am wondering - instead of looping through all features of shp1 and shp2 - if there is an easier/faster way to achieve this?

from osgeo import gdal
import ogr

vector1 = ogr.Open(path_to_shp1)
vector2 = ogr.Open(path_to_shp2)
layer1 = vector1.GetLayer()

for i in range(0, layer1.GetFeatureCount()):
    feature1 = layer1.GetFeature(i)
    if str(feature1.name) == 'MyFeatureName':
        vectorGeometry1 = feature1.GetGeometryRef()
        layer2 = vector2.GetLayer()
        for j in range(0, layer2.GetFeatureCount()):
            feature2 = layer2.GetFeature(j)
            vectorGeometry2 = feature2.GetGeometryRef()
            intersection = vectorGeometry1.Intersection(vectorGeometry2)
            wkt = intersection.ExportToWkt()
            if str(wkt) != "GEOMETRYCOLLECTION EMPTY" and str(wkt) != "POLYGON EMPTY":
                print(feature2.name)

Answer 1

I do not think that you will avoid a nested loop. A simpler way might be:

1. Open the shapefiles and copy the feature id and the geometry in two dictionaries as (id, geometry) key-value.
2. Close the shapefiles
3. Iterate over the key,value of the reference dictionary and check if the value of it intersects one of the values of the other dictionary. If no, remove this item from the dictionary.
4. You will end with the first dictionary containing all the ids of intersecting objects

Alternatively, in your solution instead of printing the features I would add an int column which will indicate if the object has 0 or more intersections.

Hope it gives you another view of the solution.