2

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.

  • I think it should be possible doing that using PyQGIS, buy you need to specify what do you mean with "location": a pair of coordinates (which format?), a cell of a (squared, hexagonal,...) vector grid, or what else? – mgri Feb 10 '17 at 13:23
  • I have edited the question and specified what I mean by location. So your advice that the python console and its functionality would be the best way to go about in doing this? Is that your answer? Why? – J.A.Cado Feb 10 '17 at 16:45
  • I suggest no longer using the "best" term because your question could be flagged as "primarily opinion-based". I said that a solution could be found using PyQGIS, but I don't know if it will be the "best" solution for you. So, the input would be a pair of coordinates? – mgri Feb 10 '17 at 16:54
  • Good point, words such as best are best avoided ;) but indeed a pair of coordinates. However , I am interested in opinions, I don't expect there to be a perfect answer but I am interested in why you would choose one specific option. – J.A.Cado Feb 10 '17 at 17:10
  • I think it would be enough using PYQGIS, especially if you want to do the comparison by using a pair of coordinates. I'll try to post my solution as an answer as soon as possible. – mgri Feb 10 '17 at 17:15
2

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:

enter image description here

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

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.