Thursday, 4 October 2018

pyqgis - QgsfeatureRequest - Store result indexes


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):


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


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...