6

I would like to find on which side of the line (left or right) my points are located with respect to the direction of the line (like knowing the rivers bank along a meandering river). But, if there are several line, I am only interested in the closest line.

I know that I could use one sided buffers like it is suggested in this post, but this doesn't work well when several lines are close to each other near the point.

For instance, there are sometimes two (or more) lines close to my point of interest, which then create two (or more) overlapping polygons covering this point. In this case, the point could be on the left of one line and on the right of the other line.

I only want to know the side of the closest road to the point (avoiding multiple buffers and selections because the database is very large. A single query maybe)

  • How about this approach gis.stackexchange.com/questions/156578/…? – user30184 Sep 7 '18 at 10:18
  • Maybe I should try PostGIS, but I was looking how to do it in QGIS. Note that I know how to do it in ArcGIS (positive or negative distance with locatePointAlongRoute) but I am looking for a solution within QGIS framework without importing my dataset into PostGIS. – radouxju Sep 7 '18 at 11:18
2

Not a single query but the following code can be executed in the Python Console. The code essentially does the following:

  • Iterates each line feature and calculates the distance from this to a selected point feature;
  • Stores these in a dictionary;
  • A proximity is used to limit iterating through each line feature (i.e. a value of 100 means only lines within 100m of the point will be included);
  • Sorts the dictionary and selects the line feature with the minimum distance from the point;
  • Gets the coordinates of the start and end vertices;
  • Determines if the point lies on the left- or right-hand side of the line (assuming the line is a single straight line).

pointLayer = QgsProject.instance().mapLayersByName('point')[0]
lineLayer = QgsProject.instance().mapLayersByName('line')[0]

# Set proximity limit
proximity = 1
distance_dict = {}

for point in pointLayer.selectedFeatures():
    start_point = point.geometry().asPoint().x()
    end_point = point.geometry().asPoint().y()
    for line in lineLayer.getFeatures():
        distance = point.geometry().distance(line.geometry())
        if distance < proximity:
            distance_dict[line.id()] = distance
    closestFeature = sorted(distance_dict.keys(), key=(lambda key:  distance_dict[key]))
    lineLayer.selectByIds([closestFeature[0]])
    for line in lineLayer.selectedFeatures():
        start_geom = line.geometry().asPolyline()[0]
        end_geom = line.geometry().asPolyline()[1]
    d = ((start_point - start_geom.x()) * (end_geom.y() - start_geom.y())) - (((end_point - start_geom.y()) * (end_geom.x() - start_geom.x())))
    if d < 0: 
        print('Left')
    elif d > 0:
        print('Right')
    else:
        print('On line')

Example of selected points (highlighted in yellow):

  • Right of line -

    Right of line

  • Left of line -

    Left of line

  • Where is the "left or right" part? – user30184 Sep 7 '18 at 11:57
  • @user30184 - Good point, forgot to add that logic in. Will edit my post soon thanks :) – Joseph Sep 7 '18 at 12:44
  • I first tought that it was working, but in fact not in the case of polyline, because it only take the first segment – radouxju Sep 10 '18 at 8:55
  • Also, a small change is needed concerning the loop: it is better to reinitialize the dictionnary for each point, otherwise there is a small risk that "far away" (> proximity) lines get used with some points. – radouxju Sep 10 '18 at 9:12
  • @radouxju - Thanks for your comment, a polyline situation would be more complex. I'll edit this post if I can come up with something :) – Joseph Sep 11 '18 at 9:15
1

@Joseph's answer is great for segment and has implemented an algorithm. However, I've found that clostsSegmentWithContext() could be used to know if the point is on the left (result[3]=-1) or on the right (result[3]=1). Note that I first screen out the possible closest lines with spatial index, which returns the closts lines based on boxes. Therefore I a a loop on the 4 closest lines, more lines could be necessary in more complex situations (but of cours it would take longer). If there are many point but only a few lines, the loop can be done systematically on each line (without "preselecting" with spatial index)

pointlayer = QgsProject.instance().mapLayersByName('p')[0] #iface.mapCanvas().layer(0)
lineLayer = QgsProject.instance().mapLayersByName('lines')[0] # iface.mapCanvas().layer(1)

testIndex = QgsSpatialIndex(lineLayer)

pointlayer.startEditing() #edit is necessary to be able to update the points

feats=[]

for p in pointlayer.getFeatures():
    mindist = 10000. #initialize a large mindist value
    near_ids = testIndex.nearestNeighbor(p.geometry().asPoint(),4)
    features = lineLayer.getFeatures(QgsFeatureRequest().setFilterFids(near_ids))
    for tline in features: #loop on 4 closest line candidtates
        closeSegResult = tline.geometry().closestSegmentWithContext(p.geometry().asPoint())
        if mindist > closeSegResult[0]:
            mindist = closeSegResult[0]
            p['leftright'] = closeSegResult[3] #update the side of line in 'leftRight' field for the closest segment 
    pointlayer.updateFeature(p)

pointlayer.commitChanges() #close edit session to flush changes

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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