[GIS] Select polygon features that intersect a line feature using pyqgis

intersectionoverlapping-featurespyqgis

How do I select all features in a polygon layer that intersect (overlap) a line feature?

I need to implement this in python using QGIS, thus pyqgis. This code illustrates what I am looking for (line 10 is sort of pseudo-code, i.e. what I am looking for):

 1 area_layer_file = "./Areas.shp"
 2 line_layer_file = "./Lines.shp"
 3 
 4 area_layer = QgsVectorLayer(area_layer_file, "Areas", "ogr")
 5 line_layer = QgsVectorLayer(line_layer_file, "Lines", "ogr")
 6 
 7 areas = []
 8 for line_feature in line_layer.getFeatures():
 9     request = QgsFeatureRequest(QgsExpression(
10             area_feature.geometry().intersects(line_feature.geometry())
11         ))
12     areas = [area for area in area_layer.getFeatures(request)]

This image illustrates my question: the red line intersects the yellow selected areas.

enter image description here

Best Answer

areas = []
for line_feature in line_layer.getFeatures():
    cands = area_layer.getFeatures(QgsFeatureRequest().setFilterRect(line_feature.geometry().boundingBox()))
    for area_feature in cands:
        if line_feature.geometry().intersects(area_feature.geometry()):
            areas.append(area_feature.id())

area_layer.select(areas)

Related Solutions

Shapefiles – Intersect Using Shapely in Python

While I'm a big user of both shapely and fiona, I wouldn't go this approach. This is a task of writing an effective SQL statement.

Using ogr2ogr with an SQLITE dialect, you can process this from a command line. Change directory to one before the shapefiles, so that all of the shapefiles are in one directory called data. OGR treats directories of shapefiles as datasets. Now try this command:

ogr2ogr -sql "SELECT ST_Intersection(A.geometry, B.geometry) AS geometry, A.*, B.* FROM nyct2010 A, h1_buf B WHERE ST_Intersects(A.geometry, B.geometry)" -dialect SQLITE data data -nln h1_buf_int_ct2

This makes h1_buf_int_ct2.shp with the same attributes as A and B, showing the intersections. You could also have done the buffering operation somewhere in there too.

The same SQL statement works with PostGIS too, so you know.

Here's how to do the same from within Python using ogr:

from osgeo import ogr
ogr.UseExceptions()
ogr_ds = ogr.Open('/path/to/data', True)  # Windows: r'C:\path\to\data'
SQL = """\
    SELECT ST_Intersection(A.geometry, B.geometry) AS geometry, A.*, B.*
    FROM nyct2010 A, h1_buf B
    WHERE ST_Intersects(A.geometry, B.geometry);
"""
layer = ogr_ds.ExecuteSQL(SQL, dialect='SQLITE')
# copy result back to datasource as a new shapefile
layer2 = ogr_ds.CopyLayer(layer, 'h1_buf_int_ct3')
# save, close
layer = layer2 = ogr_ds = None

PyQGIS – Adding and Filtering Layer with Map Canvas Extent

You can create a temporary layer based on your Postgis layer.

You can either filter features of the layer wich intersect the canvas extent, or filter the features wich are completly within the canvas extent.

On the example below, I have a layer (red), a canvas extent (red line) zoomed on the green feature, features intersecting the extent (pink features) and features within the extent (purple features).

To create those 2 layers I used this code :

import tempfile
tf = tempfile.TemporaryDirectory()

# IMPORT YOUR LAYER

if self.ctrs_cibles.featureCount() > 0:
    # Creation of the layer with the features wich intersect the extent
    self.ctrs_cibles.selectByRect(extent)
    writer = QgsVectorFileWriter.writeAsVectorFormat(self.ctrs_cibles, tf.name +'/'+'filtered_layer_intersect.shp', "utf-8", self.ctrs_cibles.crs(), "ESRI Shapefile", onlySelected=True)
    iface.addVectorLayer( tf.name +'/'+'filtered_layer_intersect.shp', '', 'ogr')
    # Creation of the layer with the features wich are within the extent
    inside_extent_ids = []
    for feature in self.ctrs_cibles.getSelectedFeatures():
        if feature.geometry().within(QgsGeometry().fromRect(extent)):
            inside_extent_ids.append(feature.id())
    self.ctrs_cibles.selectByIds(inside_extent_ids)
    writer = QgsVectorFileWriter.writeAsVectorFormat(self.ctrs_cibles, tf.name +'/'+'filtered_layer_within.shp', "utf-8", self.ctrs_cibles.crs(), "ESRI Shapefile", onlySelected=True)
    iface.addVectorLayer( tf.name +'/'+'filtered_layer_within.shp', '', 'ogr')

Best Answer

Related Solutions

Shapefiles – Intersect Using Shapely in Python

PyQGIS – Adding and Filtering Layer with Map Canvas Extent

Related Question