5

I have a cluster of points as seen below from a shapefile layer.

enter image description here

I need an automated procedure (in PyQGIS or any) of filtering out those points that are outside or further-away from the center of the clustered points, by creating a new attribute field for the shapefile layer.

The output would be the same shapefile but with an added field (named "Outlier"), where all the outliers points will be marked with a "1" and the rest with a "0".

That is the points in red circle will have value of "0" while those outside the red mark will have the value of "1".

Code Edit:

from PyQt4.QtCore import QVariant

point_layer = iface.activeLayer()
data_provider = point_layer.dataProvider()

count_ft = data_provider.featureCount()

# Create Outlier attribute field
data_provider.addAttributes([QgsField("Outliers",  QVariant.Int)])
point_layer.updateFields()

# Get points average coordinates
features = point_layer.getFeatures()
for f in features:
    #calculate average xy coordinate
    pass

# Buffer from average xy coordinate to get outliers

# update "outlier" attribute field
## if point in buffer, update with 0 else update 1
mgri
  • 16,159
  • 6
  • 47
  • 80
Umar Yusuf
  • 748
  • 1
  • 5
  • 15
  • What does your code so far look like? If you have none these Q&As may help you make a start: https://gis.stackexchange.com/search?q=%5Bpyqgis%5D+cluster – PolyGeo Oct 18 '17 at 05:32

1 Answers1

3

You may follow this procedure:

  1. Choose a proper radius value and then run the "Fixed distance buffer" algorithm (leave the "Dissolve result" option checked);
  2. Run the "Multipart to singleparts" algorithm;
  3. Iterate over each polygon features and find the outliers and, if it contains only one (or a different value) point, then it is an outlier.
  4. If the current point is an outlier, assign 1, otherwise assign 0.

The above procedure is implemented in the following script:

##Points=vector point
##Radius=number 5
##Treshold=number 1

from qgis.core import *
from PyQt4.QtCore import QVariant
import processing    

point_layer = processing.getObject(Points)
data_provider = point_layer.dataProvider()

count_ft = data_provider.featureCount()

# Create Outlier attribute field
data_provider.addAttributes([QgsField('Outliers',  QVariant.Int)])
point_layer.updateFields()
field_index = point_layer.fieldNameIndex('Outliers')

all_feats = {}
index = QgsSpatialIndex()
for ft in point_layer.getFeatures():
    index.insertFeature(ft)
    all_feats[ft.id()] = ft

buff = processing.runalg('qgis:fixeddistancebuffer', point_layer, Radius, 5, True, None)
buffered = processing.getObject(buff['OUTPUT'])

multi_to_single = processing.runalg('qgis:multiparttosingleparts', buffered, None)
mts = processing.getObject(multi_to_single['OUTPUT'])

point_layer.startEditing()
for feat in mts.getFeatures():
    geom = feat.geometry()
    idsList = index.intersects(geom.boundingBox())
    for id in idsList:
        tmp_geom = all_feats[id].geometry()
        if not tmp_geom.intersects(geom):
            del idsList[idsList.index(id)]
    if len(idsList) <= Treshold:
        for id in idsList:
            point_layer.changeAttributeValue(id, field_index, 1)
    else:
        for id in idsList:
             point_layer.changeAttributeValue(id, field_index, 0)
point_layer.commitChanges()

The script requires three inputs:

  • Points is the point vector layer you want to use;
  • Radius is the radius of the buffer;
  • Threshold is the threshold value to use for identifying the outliers (you should leave it set to 1, unless you need to manage particular situations).
mgri
  • 16,159
  • 6
  • 47
  • 80
  • I like using ##Points=vector point which populates the combobox with only point layers :) – Joseph Oct 20 '17 at 14:53
  • 1
    @Joseph me too, but I often forget to add that word! :) I have just added it, thanks. – mgri Oct 20 '17 at 14:56