Determining on which side of line points are located using QGIS

Question

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 '19 at 17:11
Hi CSK, I saw that post earlier but I could not reproduce that code on Python. — Marcel Campion, Jan 15 '19 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 '19 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 : https://gis.stackexchange.com/questions/308440/clustering-30m-pixels-into-120m-pixels-fishnet-using-qgis — Marcel Campion, Jan 15 '19 at 17:42
You can find also other answers for similar questions https://gis.stackexchange.com/search?q=line+which+side — user30184, Jan 15 '19 at 18:01
If the line is straight and has no vertices in between the start and end point you could use a trigonometric function to determine where the point falls in relation to the start and end point. Assuming you are using a Cartesian coordinate system. — Pete, Jan 15 '19 at 20:57
@Marcel Campion 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. Please, see my Editing Note. — xunilk, Jan 18 '19 at 00:40

score 14 · Accepted Answer · edited Nov 26 '23 at 17:21

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 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].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 &gt; 0:
    pos = 0
if az &lt; 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 &gt; 0:
    pos = 0
    ids_0.append(i)

if az &lt; 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! — Marcel Campion, Jan 17 '19 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: https://nathanw.net/2013/01/04/using-a-qgis-spatial-index-to-speed-up-your-code/ and try to adapt it. — xunilk, Jan 17 '19 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() — Marcel Campion, Jan 17 '19 at 14:54
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. — xunilk, Jan 17 '19 at 21:34
Hi xunilk, thank you again for your help. I am not able to use indexes lists to produce memory layer with QgsFeatureRequest. As a second best I am thinking in splitting my original point layer into 15 or 20 smaller point layers and use the first script you wrote as an answer to my question. Thank will take some time but it should work. — Marcel Campion, Jan 18 '19 at 13:45
provides an easy way to create a memory layer with the result https://gis.stackexchange.com/questions/309166/qgsfeaturerequest-store-result-indexes — Marcel Campion, Jan 18 '19 at 23:15

Babel · Answer 2 · 2024-01-11T22:46:20.320

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:

Taras · Answer 3 · 2023-12-03T18:26:35.043

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
provider.addAttributes([QgsField("position", QVariant.String)])
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:

Difference between QgsPoint, QgsPointXY and QgsGeometry.fromPointXY() in PyQGIS

score 0 · Answer 4 · answered Feb 27 '24 at 20:58

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 )

Determining on which side of line points are located using QGIS

4 Answers4

Variant

referring to the original point and line layers

accessing point layer provider

creating a field for the output

getting one single line feature

extracting line first and last point

looping over each point feature

Linked