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)
Answer
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):
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.
No comments:
Post a Comment