# Determining on which side of line points are located using QGIS

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 a value of 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 values for points located on the right of the line.

• Have you seen this question? Perhaps you can convert the SQL coded answer to Python.
– csk
Jan 15, 2019 at 17:11
• Hi CSK, I saw that post earlier but I could not reproduce that code on Python. Jan 15, 2019 at 17:32
• One of the rules on this site is that when you ask for help with coding, you must at least attempt to write some code. This is not a free coding service where people will just write code for you from scratch. So I recommend you edit your question to include a sample of the code that you've attempted so far. Otherwise your question will likely remain unanswered.
– csk
Jan 15, 2019 at 17:36
• Thank for this helpful comment. If you look at the other questions I have asked on this website you will see that I provide the code I write and I answer my own questions when I find a solution : gis.stackexchange.com/questions/308440/… Jan 15, 2019 at 17:42
• You can find also other answers for similar questions gis.stackexchange.com/search?q=line+which+side Jan 15, 2019 at 18:01

By using PyQGIS 2 (I supposed for your tag ) you can do that with `closestSegmentWithContext()` method from `QgsGeometry` and `azimuth()` method from `QgsPoint`.

I used the shapefile of the 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].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 the Python Console of QGIS, it was added the `"position"` field with a value of 0 for left and 1 for right as expected.

Editing Note:

The Process which requires more time is writing position values to the 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 list when the condition is met. The 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 the above code, I generated 900,000 random points and ran it. The result was produced in less than a second and, on average, each list has about 450,000 points as expected (see the below image).

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

• xunilk, this solution is very good. My point layer is huge (950.000 obs) so it takes a long time to run but it actually worked so thank you! Jan 17, 2019 at 13:06
• Most welcome, glad it helped. So, I didn't know how many points had your point layer. If you want to speed up your code you can use a QGIS spatial index. Please, see this post: nathanw.net/2013/01/04/… and try to adapt it. Jan 17, 2019 at 13:54
• Hi again xunilk, I am not sure about how to speed up your code. I am not proficient in pyqgis so it's like making sentences in a foreign language. Here is what I have done so far: def withindex(): index = QgsSpatialIndex() for feat in points_feats: index.insertFeature(f) for feature in allfeatures.values(): ids = index.intersects(feature.geometry().boundingBox()) for id in ids: if f==feature:continue pt = feat.geometry().asPoint() Jan 17, 2019 at 14:54
• Ok. I will see later. Jan 17, 2019 at 18:06
• I spent some time to add one Spatial Index but runnig time is similar in modified code, So, I found out that this approach is useful when there are two loops where you check each feature against every other feature, In your code there is only one loop. Approach must be different. Jan 17, 2019 at 21:34

You can also use QGIS expressions for this. See below for the expression to get a 0 or 1 output for each point, depending on which side of the line it is (see below for curved lines). The solution is based on function `azimuth()` to calculate the north-based angle of the line connecting the closest point on the line to the current point:

``````with_variable(
'line',  -- name of your line layer
overlay_nearest ('line',\$geometry)[0],
azimuth (start_point(@line), end_point(@line)) >
azimuth(
closest_point (@line,\$geometry),
\$geometry
)
)
``````

Point layer labeled with the expression above:

This works for a straight line. If you have a curved line ("on which side of the rivers are my points of interest?"), use this expression instead:

``````with_variable(
'ln',
overlay_nearest ('line',\$geometry)[0], -- name of your line layer
with_variable (
'pa',
degrees (azimuth(closest_point (@ln,\$geometry),\$geometry)),
with_variable (
'la',
line_interpolate_angle(
@ln,
line_locate_point (@ln,closest_point (@ln,\$geometry))
),
case
when (@la < 90 and @pa > 270) then 0
when (@la > 270 and @pa <90 )then 1
else @la < @pa
end
)))
``````

Points labeled with the expression and a variant of it to assign different colors:

#### Variant

As @Taras asked in a comment: what happenes if the point is exactly on the line? You can include a condition (e.g.return `on line`) for that case by modifying the expression with an additional if-condition:

Version for straigth line:

``````with_variable(
'line',  -- name of your line layer
overlay_nearest ('line_straigth',\$geometry)[0],
if(
intersects (buffer (\$geometry,0.0001), @line),
'on line',   -- result when point is on the line
azimuth (start_point(@line), end_point(@line)) >
azimuth(
closest_point (@line,\$geometry),
\$geometry
)
)
)
``````

Version for curved line:

``````with_variable(
'ln',
overlay_nearest ('line',\$geometry)[0], -- name of your line layer
with_variable (
'pa',
degrees (azimuth(closest_point (@ln,\$geometry),\$geometry)),
with_variable (
'la',
line_interpolate_angle(
@ln,
line_locate_point (@ln,closest_point (@ln,\$geometry))
),
if (
intersects (buffer (\$geometry,0.0001), @ln),
'on line',  -- result when point is on the line
case
when (@la < 90 and @pa > 270) then 0
when (@la > 270 and @pa <90 )then 1
else @la < @pa
end
))))
``````

Result with an additional few points exactly on the line:

In PyQGIS 3 there is a handy method available, namely `segmentSide()` from the `QgsGeometryUtils` class.

Let's assume there are two layers: a line layer with a single feature, and a point layer with nine features in it.

I am working with the `EPSG:25832`.

Proceed with `Plugins > Python Console > Show Editor` and paste the script below:

``````# imports
from qgis.core import QgsProject, QgsField, QgsPoint, QgsGeometryUtils
from PyQt5.QtCore import QVariant

def find_position(start_point: QgsPoint, end_point: QgsPoint, test_point: QgsPoint) -> str:
"""
Finds out on which side of a line a point is located
Parameters:
==========
:param start_point: the line's starting point
:param end_point: the line's ending point
:param test_point: the testing point
"""
positions = {'left': -1, 'right': 1, 'on the line': 0}
position = QgsGeometryUtils.segmentSide(start_point, end_point, test_point)
result = list(positions.keys())[list(positions.values()).index(position)]
return result

# referring to the original point and line layers
point_layer = QgsProject.instance().mapLayersByName("deg_centroids")[0]
line_layer = QgsProject.instance().mapLayersByName("deg_line")[0]

# accessing point layer provider
provider = point_layer.dataProvider()
# creating a field for the output
point_layer.updateFields()

# getting one single line feature
line = line_layer.getFeature(0).geometry()
# extracting line first and last point
first_vertex = line.vertexAt(0)
last_vertex = line.vertexAt(1)

point_layer.startEditing()

# looping over each point feature
for feature in point_layer.getFeatures():
# getting each feature geometry as QgsPoint
point = QgsPoint(feature.geometry().asPoint())
# providing new values for the "position" attribute
feature["position"] = find_position(first_vertex, last_vertex, point)
point_layer.updateFeature(feature)

point_layer.commitChanges()
``````

Change the inputs in the `point_layer` and `line_layer` variables. Press `Run script` and get the output that will look like this:

References:

For the case of a roadway where we used "CL", "LT", or "RT" to specify the location of a point, the last part of the statement by Babel could be modified as below: if (intersects (buffer (\$geometry,0.0001), @ln), 'CL', -- result when point is on the line case when (@la < 90 and @pa > 270) then 'LT' when (@la > 270 and @pa <90 )then 'RT' when (@la > @pa) then 'LT' else 'RT' end )