5

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.

enter image description here

7
  • 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
  • 1
    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
  • 1
    You can find also other answers for similar questions gis.stackexchange.com/search?q=line+which+side
    – user30184
    Jan 15, 2019 at 18:01

4 Answers 4

14

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:

enter image description here

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

enter image description here

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

enter image description here

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.

7
  • 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
  • 1
    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.
    – xunilk
    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
  • 1
    Ok. I will see later.
    – xunilk
    Jan 17, 2019 at 18:06
  • 1
    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, 2019 at 21:34
8

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: enter image description here

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: enter image description here


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: enter image description here

0
7

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.

input

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 run script and get the output that will look like this:

result


References:

0

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 )

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.