qgis - Determine on which side of a line points are located

I have a point shapefile (point.shp) and a line shapefile (line.shp). I want to determine on which side of the line each point is located.

The best would be to use a pyqgis command to add a field to the point.shp with value 0 for left and 1 for right.

An alternative that would work is to determine the nearest distance between each point and the line with distance taking negative values for points located on the left and positive value for points located on the right of the line.

Answer

By using PyQGIS 2 (I supposed for your tag pyqgis) you can do that with 'closestSegmentWithContext' method from QgsGeometry and 'azimuth' method from QgsPoint. I used shapefile of following image:

to test below code:

from PyQt4.QtCore import QVariant


registry = QgsMapLayerRegistry.instance()

points = registry. mapLayersByName('points_grid')
line = registry. mapLayersByName('line')

position_field = QgsField( 'position', QVariant.Int)

points[0].dataProvider().addAttributes([position_field])
points[0].updateFields()


points_feats  = [ feat for feat in points[0].getFeatures() ]
line_feat = line[0].getFeatures().next()

idx = points[0].fieldNameIndex('position')

points[0].startEditing()

for feat in points_feats:
    pt = feat.geometry().asPoint()

    point_in_line = line_feat.geometry().closestSegmentWithContext(pt)[1]
    az = pt.azimuth(point_in_line)

    if az > 0:
        pos = 0
    if az < 0:
        pos = 1

    feat[idx] = pos
    points[0].updateFeature( feat )


points[0].commitChanges()

After running it at Python Console of QGIS, it was added 'position' field with value 0 for left and 1 for right as expected.

Editing Note:

The Process which requires more time is writing position values to respective field. To avoid that is preferable to use an enumerate list for points features and print point indexes (for left and right points) into a separate lists when the condition is met. Following code can do that.

from PyQt4.QtCore import QVariant

registry = QgsMapLayerRegistry.instance()


points = registry. mapLayersByName('random_points_900000')
line = registry. mapLayersByName('line')

points_feats  = [ feat for feat in points[0].getFeatures() ]
line_feat = line[0].getFeatures().next()

idx = points[0].fieldNameIndex('position')

ids_0 = []

ids_1 = []

for i, feat in enumerate(points_feats):
    pt = feat.geometry().asPoint()
    point_in_line = line_feat.geometry().closestSegmentWithContext(pt)[1]
    az = pt.azimuth(point_in_line)

    if az > 0:
        pos = 0
        ids_0.append(i)


    if az < 0:
        pos = 1
        ids_1.append(i)

print "Wait..."
print "Time: %s seconds " % timeit.timeit()

To test above code, I generate 900,000 random points and ran it. Result was produced in less than a second and, on average, each list has about 450,000 points as expected (see below image).

You can use indexes lists for producing memory layers with a QgsFeatureRequest or introduce them in point layer by using a more efficient method.