1

using QGIS 2.18.16 I am trying to figure out how to find out what polygons in layer one (blue) intersects the polygons in layer 2 (red\pink). I know that there will be more than one intersection for most of the polygons in layer 1 (blue) but not sure how to calc it. In a perfect world I would like to add it as a attribute to the polygon in layer 1.

enter image description here

  • what data system are you using? – DPSSpatial May 30 '18 at 19:31
  • Data System? I am using QGIS 2.18 and shape flies – lowsparked May 30 '18 at 19:49
  • I have used the intersection function under Vector> Geoprocessing Tools . This gives me the data but not in one to many way i am looking for it. it is many to one. Any help would be great. – lowsparked Jun 1 '18 at 17:26
  • I would prefer to do this in a relational databases where many-to-one relationships are better handled, especially for geometry... you may be able to use the DB manager in QGIS to write a SQL spatial query against your shapefiles to join the features - if that might work for you, I'll see if I can write a sample... – DPSSpatial Jun 1 '18 at 17:45
  • DISCLAIMER: not a QGIS guy, but GIS is GIS, right? When you use the Intersection Function under Vector>Geoprocessing Tools, switch the source and target parameters so you are intersecting red 1/16 section polys with blue plat maps ... I mean reverse the two arguments of the function. Think B on A, not A on B. -- This isn't an "Answer" to your problem but rather assistance with the concept behind the QGIS function. – JasonInVegas Jun 14 '18 at 22:52
2
+25

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):

quarter sections and layer 1

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".

settings for join by location

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

duplicated polygons

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()

results after python code

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.

2

I´d go for a Virtual Layer using the natively supported SQLite/SpatiaLite SQL dialect:

  • Layer | Add Layer | Add/Edit Virtual Layer...
  • Embedded Layers | Add (from file) / Import (from project) | select both polygon layers

  • Assuming their names are layer1 and layer2 with an id column each and an attribute val in layer2 that you are interested in, running

    SELECT a.*,
           GROUP_CONCAT(b.val) AS vals
    FROM layer1 AS a
    JOIN layer2 AS b
      ON ST_Intersects(a.geometry, b.geometry)
    GROUP BY a.id
    

    in it´s Query field will copy all attributes/geometries from layer1 plus a concatenated string field vals, with the val's from layer2 that intersect each polygon in layer1, into the Virtual Layer.

  • Click Add and save to file (e.g. layer context | Save as...)



The more versatile way would be to use the DB Manager (as @DPSSpatial also mentioned in his comment I only just read...); this is a full scale DB client, capable of also runnning queries against the project layers (you'll find all currently loaded layers under Virtual Layers | QGIS Layers on the left) using the SQLite/SpatiaLite SQL dialect;

  • E.g. when executing the same query, you will be presented a result table including a geometry field to check against.
  • Tick Load as new layer and save the created layer to file (don´t forget to specify a geometry column if not already set).
2

If you are not familiar with or just not willing to use SQL in any way, you can also simply use the Field Calculator (the summary aggregate functions should be available from QGIS 2.16.x on):

  • layer1 | Open Field Calculator | New Field
    • Output field type = Text (string)
  • type in Expression:

    aggregate(layer:=<layer>, aggregate:='concatenate', concatenator:=',', expression:="<field>", filter:=intersects($geometry, geometry(@parent)))
    
    • select layer2 name as <layer> under Map layers
    • replace <field> with fieldname of interest in layer2 (double quotes have to stay!)
    • if <field> is not of type string, use expression:=to_string("<field>")
  • OK should add the above field to the attribute table of layer1 as a concatenated string

  • save/quit edit mode (has been auto enabled)

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.