1

I recently posted a question. My original goal was to find a way to determine on which side of a line points are located. @xulnik answered this question with two codes. One returns a binary [0;1] column indicating on which side points are located. However as the number of points for which I want to know whether they are on the right or the left of a given line is huge (about 1,490,000) @xulnik proposed an alternative based on enumerate list for points features point indexes.

In his answer @xulnik proposes that I use indexes lists to produce a memory layer with QgsFeatureRequest. I have tried things based on @gsherman 's answer to this post but did not succeed. Basically I am trying to find a way to either export the result of the following code to a csv file or to create a layer that contains the result of the code.

from PyQt4.QtCore import QVariant

registry = QgsMapLayerRegistry.instance()

points = registry. mapLayersByName('pixels_id')
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)
2

For producing a memory layer with a request by using left points (ids_0) from line, you have to include several lines of code (see #added lines mark below):

from PyQt4.QtCore import QVariant
import timeit

registry = QgsMapLayerRegistry.instance()

points = registry. mapLayersByName('random_points_10000')
line = registry. mapLayersByName('line')

points_feats  = [ feat for feat in points[0].getFeatures() ]
line_feat = line[0].getFeatures().next()

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)

#added lines

epsg = line[0].crs().postgisSrid()

uri = "Point?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'request',
                           'memory')

prov = mem_layer.dataProvider()

request = QgsFeatureRequest().setFilterFids(ids_0)

iter = points[0].getFeatures(request)

feats = [ feat for feat in iter ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

print "Wait..."
print "Time: %s seconds " % timeit.timeit()

I tried above code out with shapefiles of following image (points layer has 10,000 random points):

enter image description here

After running above code at Python Console of QGIS, it was produced a memory layer with left points (ids_0) from line as it showed below. It worked.

enter image description here

  • this code works perfectly (and even with a huge number of points). Great! – Marcel Campion Jan 18 at 23:14
  • You're welcome. – xunilk Jan 18 at 23:24

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.