2

I have 2 polygon shapefiles: A (with integer field "a_id") and B (integer field "b_id").

My goal is to extract intersections of polygons in layer A with polygons in layer B where attribute "a_id" is equal to attribute "b_id" + 55.

I'm trying to implement it with the QGIS graphical modeler, but "Select by expression" algorithm supports only one layer for running expressions, and "Extract by location" has no attribute filter.

2
  • Could you compute all intersections first and then check the attributes?
    – underdark
    Commented Aug 1, 2016 at 19:54
  • @underdark I've tried it but calculating all intersections is quite computation-intense on a large volumes of data and it felt like there should definitely exist easier ways to do the task.
    – Basile
    Commented Aug 3, 2016 at 10:13

1 Answer 1

2

The best solution for me was to open directory containing both layers as a single GDALDataset object, make an OGR SQL query to it and save the result to a new layer. It allows doing spatial and attribute queries simultaneously and can be done with the QGIS console:

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.a_id, B.b_id
    FROM A, B
    WHERE ST_Intersects(A.geometry, B.geometry)
    AND CAST(A.a_id as INTEGER) = CAST(B.b_id as INTEGER) + 55;"""
layer = ogr_ds.ExecuteSQL(SQL, dialect='SQLITE')
# copy result back to datasource as a new shapefile
layer2 = ogr_ds.CopyLayer(layer, 'result_shapefile')
# save, close
layer = layer2 = ogr_ds = None

The idea was picked here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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