Optimal Site Selection Vector Files

Question

I have 5 different polygon vectors:

• All have over 20.000 features;
• All have multiple attributes;
• All overlap each other at least more than 80%, but are shaped differently;
• All display information of large areas, about 1/3 of the North Sea;
• All have an attribute called "benefit" which ranges between 0-1 (I normalized the data);

My mission: Using the attribute "benefit" I want to find optimal locations for a specific activity. The sum of attributes for 5 different vector files, would indicate where I can find such optimal sites. The higher the better.

For example, let's say I choose two different geographic locations close to each other somewhere off the eastern coast of England. For each location I would retrieve the relevant value of the attribute "benefit".

For the attribute "benefit" I would get:

Location 1
layer 1: 0.5
layer 2: 0.3
layer 3: 0.6
layer 4: 0.8
layer 5: 0.9 

Location 2
layer 1: 0.1
layer 2: 0.3
layer 3: 0.2
layer 4: 0.1
layer 5: 0.4

Location 1 would be better compared to location 2. 
Sum layer 1: 3.1 vs Sum layer 2: 1.1

THE QUESTION: What would be the best methdos to do this for all layers across the all extends?

At my disposal I have: QGIS, Python and pgAdmin

1) QGIS would be the least technical but I'm not even sure it's possible. With such large vector files I think it would crash. I could use the python console maybe...?

2) pgAdmin might be quicker, but for each feature I would need to check if it overlaps any of the other features...and with 5 vector files each containing more than 20,000 features...(some even 100,000)...I'm not sure if this is possible.

3) Python might also allow for a quick analyses of the data. But I'm not sure which packages to use.

Answer 1

Context

On the basis of your request, I initially thought about a solution using the Python Console: starting from the inserting of a pair of coordinates, you would be able to get a print in the Console containing the values stored on each intersected polygon and finally the total sum for that location, as in your example.

This would be a simple and logic way to approach the problem, but not really efficient when dealing with more than a couple of locations (for several programming reasons, but this is not the place for discussing them).

In my opinion, a more performing solution would be the using of a point layer as input: in this way, you may literally increase the speed of your analysis and avoiding the manual insertion of the coordinates.

Solution

I propose a solution using a PYQGIS script. You only need one "ingredient" before starting the analysis, i.e. a point layer containing your coordinates: this is a very simple preliminary task and there are several ways for doing this, but I can't explain it here because my answer would become off-topic (remember also that you only need to create some points as geometries, don't worry about the attributes).

Once you have created a point layer, you may run this script from the Processing Toolbox:

##Points=vector point
##polygons=multiple vector

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

layers = polygons.split(';')
points = processing.getObject(Points)
fields = points.pendingFields()

indexes = []
feats= {}
sum = 0
for poly in layers:
    layer = processing.getObject(poly)
    l_name = layer.name()
    fields.append(QgsField('%s' %(l_name), QVariant.Double, '', 10, 3))
    index = QgsSpatialIndex()
    for feat in layer.getFeatures():
        index.insertFeature(feat)
        feats[feat.id()] = feat["benefit"]

    indexes.append(index)

# This layer is the final output
crs = points.crs().toWkt()
output = QgsVectorLayer('Point?crs='+ crs, 'results' , 'memory')
prov = output.dataProvider()
fields.append(QgsField('Sum', QVariant.Double, '', 10, 3))
prov.addAttributes(fields)
output.updateFields()

for point in points.getFeatures():
    sum=0
    attrs = point.attributes()
    PtGeom = point.geometry()
    for k in indexes:    
        idsList = k.intersects(PtGeom.boundingBox())
        if not idsList:
            value = 0
        else:
            value = feats[idsList[0]]
        attrs.append(value)
        sum += value
    attrs.append(sum)
    outGeom = QgsFeature()
    outGeom.setGeometry(PtGeom)
    outGeom.setAttributes(attrs)
    prov.addFeatures([outGeom])

QgsMapLayerRegistry.instance().addMapLayer(output)

This code will create a new point memory layer, geometrically identical to the one you inserted as input, but also containing new fields which store the value of every (eventually) intersected polygon and the final sum:

(In the above image, I ran the analysis using a point layer having 8 point features. "Layer 1" to "Layer 5" are the names of some sample layers I used).

I think that in this way it would be very easy comparing the results from different location at once!