1

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

enter image description here

  • Have you seen this question? Perhaps you can convert the SQL coded answer to Python. – csk Jan 15 at 17:11
  • Hi CSK, I saw that post earlier but I could not reproduce that code on Python. – Marcel Campion Jan 15 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 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/… – Marcel Campion Jan 15 at 17:42
  • 1
    You can find also other answers for similar questions gis.stackexchange.com/search?q=line+which+side – user30184 Jan 15 at 18:01
5

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 shapefile of 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 Python Console of QGIS, it was added 'position' field with value 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 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. 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 above code, I generate 900,000 random points and ran it. Result was produced in less than a second and, on average, each list has about 450,000 points as expected (see below image).

enter image description here

You can use indexes lists for producing memory layers with a QgsFeatureRequest or introduce them in 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 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 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 at 14:54
  • 1
    Ok. I will see later. – xunilk Jan 17 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 at 21:34

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.