6

I am using QGIS 2.18.1 Las Palmas, in Windows 10, and that's the issue I want to clarify:

I've got 2 polygon shapes. 1- One with river basins, and 2- other one with different levels of dessertification risk:

River basins

dessertificationrisk

In the attribute tables for each one, I've got names for each river basin in the first map, and in the second one, I've got 5 different levels (differenced by colors in the picture)

So, what I need is not only cut or intersect both of them. What I need is to know which level of dessertification risk is more predominant, within each river basin. In other words, I need a final shape similar as the river basins shape, but with an additional column saying the level of dessertification risk which occupies the biggest % of surface within each basin.

For example: I will zoom into one of the basins, setting the river basins without filling, so you can see the dessertification risk layer behind:

enter image description here

We see different colors within the polygon. Initially, it seems that the most predominant color is purple, or could be the light green.

What I need is to determine which one covers the biggest % of surface within this polygon, and set it as a new value for the river basins polygons.

I've tried different things but I don't find the solution.

  • Most certainly the ArcGIS workflow gis.stackexchange.com/questions/217729/… can be implemented in QGIS – FelixIP Dec 1 '16 at 21:41
  • Thanks @FelixIP ! I've tried to follow the same workflow, but in step 2, I can't find the option "summarize" in QGIS. Where can I find it? Thanks! – Alex Fernández Poulussen Dec 5 '16 at 0:05
  • There is no summarise in my solution, although I'd be surprised if such function does not exist in qgis. I don't know qgis, sorry – FelixIP Dec 5 '16 at 0:28
  • Sorry, my fault. I refer to the case below yours, done by @jbchurchill . In your case, that's exactly what I need to do, but maybe because of my lack of ability with PyQgis I don't know how to implement the workflow in QGIS. – Alex Fernández Poulussen Dec 5 '16 at 10:16
  • My solution based 100℅ on tools,no scripts – FelixIP Dec 5 '16 at 17:27
5

NOTE I edited the code because the questioner reported some issues with the results.


I have created a sample dataset to reproduce the issue and I have made some assumptions:

  • the layer of the river basins (colored with light green) stores the name of each river basin in a field called "BASIN_NAME";
  • the layer with different levels of desertification risks (colored with a color ramp of reds) stores the value of the risk in a field called "RISK_LEVEL";
  • the level of desertification risk is formatted as an integer value (but this can be easily adapted to your specific needs).

enter image description here

You may use this script:

# Layer of the river basins
##risks=vector
# Layer with different levels of desertification risk
##basins=vector

from qgis.core import *
from qgis.PyQt.QtCore import QVariant

layer1 = processing.getObject(risks)
crs = layer1.crs().toWkt()
layer2 = processing.getObject(basins)

# Create the output layer
outLayer = QgsVectorLayer('Polygon?crs='+ crs, 'basins_new' , 'memory')
prov = outLayer.dataProvider()
fields = layer2.pendingFields() # Fields from the input layer
fields.append(QgsField('PREDOMIN_RISK', QVariant.Int, '', 10, 0)) # Name for the new field in the output layer
prov.addAttributes(fields) # Add input layer fields to the outLayer
outLayer.updateFields()

# Create a dictionary and a spatial index with the features from the previous intersection
allfeatures = {}
index = QgsSpatialIndex()
for feat1 in layer1.getFeatures():
    index.insertFeature(feat1)
    allfeatures[feat1.id()] = feat1["RISK_LEVEL"]

for basin in layer2.getFeatures():
    inAttr = basin.attributes() # Input attributes
    basin_geom = basin.geometry() # Input geometry
    idsList = index.intersects(basin_geom.boundingBox())
    count = 0
    req = QgsFeatureRequest().setFilterFids(idsList)
    tmp_dict = {} # Temporary dictionary containing all the features inside the current river basin
    for elem in layer1.getFeatures(req):
        temp_geometry = elem.geometry()
        if basin_geom.intersects(temp_geometry):
            itx = basin_geom.intersection(temp_geometry)
            tmp_dict[elem.id()] = itx.area() # Calculate the area
    if len(tmp_dict) > 0:
        max_key = max(tmp_dict, key=tmp_dict.get) # Evaluate the key with the maximum value of area
        inAttr.append(allfeatures[max_key])# Add the desertification risk value from the feature referring to max_key

    outGeom = QgsFeature()
    outGeom.setAttributes(inAttr) # Output attributes
    outGeom.setGeometry(basin_geom) # Output geometry
    prov.addFeatures([outGeom]) # Output feature

# Add the layer to the Layers panel
QgsMapLayerRegistry.instance().addMapLayer(outLayer)

The code will create a new polygon layer (as memory layer), having the same fields of the river basins layer, plus one additional field (called "PREDOMIN_RISK") where the predominant level of desertification risk that intersects the current river basin is stored.

For example, zooming into one of the river basins (higlighted with a blue circle):

enter image description here

we will have two different values of desertification risk (the label shows the corresponding value from the field "RISK_LEVEL"):

enter image description here

For this river basin, the attribute table from the output layer will show:

enter image description here

  • Dear @HowToInQGIS thanks soo much! the output is exactly what I'm looking for, however, I couldn't make it work. Perhaps because of my lack of ability with scripts. I've even created a copy of my layers so they have the same name (risks and basins) with the same field names ("RISK_LEVEL" and "BASIN_NAME"). Then I opened the Python console, pasted your script, and click run. The message I get is: – Alex Fernández Poulussen Dec 5 '16 at 19:36
  • 1
    The script will work if you use the Processing framework. From the main QGIS Window, go to Processing and then click on Toolbox. From the Processing toolbox (just activated on the right side of the main QGIS window), go to Create new script menu under the Tools group in the Script algorithms block of the toolbox. Then, copy & paste my code and it will magically work! You could be also run it from directly from the Python Console, but you need to change the way with which you interact with the input files. Please let me know if you encounter any other problem. – mgri Dec 5 '16 at 20:17
  • 1
    For additional referece, see also this [link] (docs.qgis.org/2.8/en/docs/user_manual/processing/scripts.html) from the QGIS Documentation. Once you will be able to work with scripts, it will be very easy to adapt the code to your specific needs (e.g. you may modify it using your original field names). – mgri Dec 5 '16 at 20:23
  • 1
    Hmm, I can't start doing some tests because I can't reproduce your issue. The best would be if you can select the features I see in the image you posted, saving them (for layerA and layerB) and then uploading them somewehere. If the intersection seems ok for you, you may also loading only the output from the intersection, so it will be easier for me doing some test. Is it ok for you? – mgri Jan 24 '17 at 14:43
  • 1
    Let us continue this discussion in chat. – mgri Jan 24 '17 at 15:35

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.