Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It only takes a minute to sign up.
Sign up to join this community
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.
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 > 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.
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:
You can also use QGIS expressions for this, using this 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:
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:
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:
