[GIS] Lines and polygons intersection with tolerance in QGIS

intersectionlinepolygonqgistolerance

I have a layer with polygons and a layer with polylines.

I am looking for the list of polygons crossed by the polylines.

Just here it's relatively simple: intersection function.

That said, the layout of the line is not perfect, so I would like to exclude polygon for which the line only touches the polygon contour or does not fit more than 1m into the polygon.

Best Answer

I'd go for SQL (QGIS natively supports the SpatiaLite 4.3.0 dialect for any layer query) and check for the spatial relation with the negative buffered polygons:

QGIS | DB Manager
- connect to Virtual Layers (double-click) in the left panel (make sure your layers are loaded in the project)
- open the SQL Window (second icon to the left)
  - in the top panel, run (Execute)
```
SELECT a.<id>,
       group_concat(b.<id>) AS its_lines,
       a.geometry 
FROM <polygon_layer> AS a
JOIN <line_layer> AS b
  ON ST_Relate(ST_Buffer(a.geometry, -1), b.geometry, 'T********')
GROUP BY a.id;
```
  - load the produced table as a layer to the project by ticking the Load as new layer checkbox and fill in the fields

The above query will return a table with all polygon <id>s (paste the actual id field name here), their geometries and a concatenated string of line <id>s (same here) that intersect, but not touch, the negative buffered polygons of your layer.

Important:
Since you asked for meters, I put -1 in ST_Buffer. However, the function will treat the passed number in CRS units; here, this implies that you are working with a CRS whose units are in meter! If you are working in LatLon, e.g. EPSG:4326, the parameter in ST_Buffer will be treated as degrees...and thus makes no sense. In that case, a transformation into an appropriate projection is necessary; you can do that with
... ST_Buffer(ST_Transform(a.geometry, <srid_of_metric_CRS>), -1) ...

Related Solutions

[GIS] One to many intersection of polygons

I found the following to work well on my small sample, though it involves repeated iteration through the layer1 features, so it might be slow on a large dataset. Also, I used QGIS 3.0.0. I am not sure how much it differs from 2.18.16.

First, load data (pardon my crooked quarter sections):

Second, use "Vector > Data Management Tools > Join Attributes By Location" to join the two layers. I selected "contains" in addition to "intersects" under geometric predicates to cover cases like polygon 2, where it doesn't cross the border of a quarter section. Use layer 1 (blue) as the input and the quarter sections as the join layer. Use the create separate feature for each located feature option under "join type".

This will create duplicate polygons wherever a layer1 polygon intersects with more than quarter section.

Lastly, we use the new joined layer to write the correct quarter section names into the original layer1 file with a bit of python.

# add layers and add new attribute to layer1
joined = QgsVectorLayer("...some path/joined.shp", "Joined layer", "ogr")
layer1 = QgsVectorLayer("...some path/layer1.shp", "layer1", "ogr")
res = layer1.dataProvider().addAttributes([QgsField("QS_list", QVariant.String)])
joined.updateFields()
layer1.updateFields()

# iterate through the layer1 features and for each,
# iterate through the joined features and for each,
# if the layer1 and joined IDs match, add the QS name to a list
# finally sort the list and write it into the new layer1 attribute
for feature in layer1.getFeatures():
    ID = feature.id()
    QS_list = []
    for feature_joined in joined.getFeatures():
        attrs = feature.attributes()
        attrs_joined = feature_joined.attributes()
        if attrs_joined[0] == attrs[0]:
            QS_list.append(str(attrs_joined[1]))
    layer1.startEditing()
    layer1.changeAttributeValue(ID, 1, str(sorted(QS_list)))
    layer1.updateFields()
    layer1.commitChanges()

There appears to be a trick to doing this in ArcGIS also (and it looks way easier). Spatial Join’s hidden trick or how to transfer attribute values in a One to Many relationship

Here are the steps (the post on Esri Australia Technical Blog also has screen captures). I didn't test this.

In the Spatial Join tool’s dialog navigate to the Zoning field in the dataset that is being joined and Right-click on it. Choose Properties.

Set the following:

Length: set the value which will be capable of storing the appropriate number of characters

Merge Rule: use Join

Delimiter: use any delimiter you like

Set the cardinality in the Spatial Join tool to “One to One” Run the tool. The resulting dataset should display zones that each of the property boundary polygons intersect, stored as text values separated by the delimiter of your choice.

Shapely – Fast Intersection of Polygons with a Bounding Polygon

To speed up this kind of queries you have to use a spatial index. Use the STRtree class to build an rtree and exponentially speed up this operation. Here is the documentation with an example:

https://shapely.readthedocs.io/en/stable/manual.html#str-packed-r-tree

EDIT:

After seen the images, the STRtree will not be usefull at all. In this case, all the polygons intersect the bounding polygon so you need a different kind of optimization. I'll give you some ideas on how you could optimize this:

When you know a geometry must be tested against all (or almost all) the geometries in a set, you can 'prepare' this geometry to speed up testing operations. It don't supports 'overlaps' but you could use 'contains' to test the polygons that falls entirely within and 'intersects' for the boundary polygons.
Knowing that the polygons you are using are convex, instead of testing if the polygon is within the boundary polygon, you could test if the centroid is within, and if not test the intersection of the full polygon. If you don't know if they whould be convex, you could also use a representative point.
I'm thinking that, in some cases, the polygons that intersect the boundary could have the centroid within. You could avoid this by using a negative buffer of the boundary polygon.
You have a 'boundary polygon' but you could also have the boundary of this 'boundary polygon', wich is a linestring that only intersects the polygons that cross that line.
Also, you can extract the segments of the boundary polygon and use the STRtree to index it. Then iterate the polygons to see if they intersect any segment.

Maybe the most optimized way I know is to test if the centroid of the polygons are within a prepared negative buffer of the boundary polygon. If the centroid is not within, seems that you don't have polygons on the outside, but im not 100% sure so I'll give you two options:

If there are no polygons on the outside, then you know that if the centroid is not within it will intersect the boundary.
Else, you could test the intersection with the 'line boundary', or if it's not fast enougth, use the STRtree of the segments.

Best Answer

Related Solutions

[GIS] One to many intersection of polygons

Shapely – Fast Intersection of Polygons with a Bounding Polygon

Related Question