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
2

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).
  • 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

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.