1

I'm pretty new to QGIS and python in general. For my work I'm faced with a bit of a problem.

I got 100-m-segments of several streams and polygon-features of the surrounding substrate. The polygon features determine if sand, gravel or stones are the primary substrate in the stream segments. So I tried using a simple spatial join. If I do that however only the segments that are 100% inside a certain substrate are correctly joined. If a segment is overlapping with two substrates, then the bigger overlapping polygon (> 50%) should be used as the primary substrate.

Is there a plugin for QGIS (version 2.18) that can do this? Otherwise, how can I create a corresponding Python script?

I searched for hours and found similar problems only with two polygon features so I tried to modify a script but failed. (Percentage of polygon in one shapefile within polygon of another)

Here's an example of overlapping substrates:stream segments and substrates

  • I tried this: 1. I intersected the 100-m-segments with the substrate, thus generating x-amount of undersegments (determined by number of polygons crossed) 2. I calculated the percentage each of these undersegments inside the polygons. 3. now I wanted to dissolve the "undersegments" (they got the same "et_id" but with different substrates and percentages). However only one set of substrates and percentages are carried over. So if a 100-m segment is intersecting 3 different substrates only the substrate at the starting point gets carried over to the attribute table. – Atlan Aug 9 '18 at 14:08
  • Check out this question/answer. After you have intersected your streams with the substrate layer, you could use the concatenate function described by Kazuhito to aggregate the substrate types and their percentages into fields before dissolving again. – Gabriel C. Nov 30 '18 at 15:16
0

First the imports:

import gdal
import ogr

Then, load the polygons

stream = ogr.Open("stream.shp")
substrates = ogr.Open("substrate1.shp")
stream_l = stream.GetLayer()
sub_l = substrates.GetLayer()

create the output file

driver = ogr.GetDriverByName("ESRI Shapefile")
outputds = driver.CreateDataSource(outputfn)
outputlayer = outputds.CreateLayer(outputfn, geom_type=ogr.wkbPolygon)

The following loop creates a with the stream based on the largest substrate present. The area of overlap, between a stream and substrate is calculated using gdal's Intersection. It is stored in a dictionary and the sorted based on these numbers. That's the idea here. I have not tested the script, so maybe you'll have to do some edits.

for stream_f  in stream_l:
    s_geom = stream_f.GetGeometryRef() 

   #Get overlap size by using a clip
    dict = {}
    for i,sub_feat in enumerate(sub_l):
        geom = sub_feat.GetGeometryRef()
        clip = geom.Intersection(stream_geom)
        dict[clip.GetArea()]=i 

    #Sort based on size and get most overlapping polygon
    sorted(dict.items(), key=lambda s: s[0])
    polynum = dict[dict.keys[0]]
    largest_feat = sub_l.GetFeature(polynum)

    #enter the fieldname at #fn
    field_largest_feat = largest_feat.GetField(#fn)

    #add this substrate feature to the original stream feature
    defn = outputlayer.GetLayerDefn()
    feature =ogr.Feature(defn)
    feature.CreateField(ogr.FieldDefn("substrate",ogr.OFTString))
    feature.SetField(str(field_largest_feat))
    feature.SetGeometry(s_geom)

    outputlayer.CreateFeature(feature)

outputds = None

If there is a problem, please ask! Hopefully this helps :)

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.