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)
Marcel Campion
  • 503
  • 3
  • 13

1 Answers1

3

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

xunilk
  • 29,891
  • 4
  • 41
  • 80