2

I am working with pyqgis in QGIS3 and i want to create a point layer from an intersection beetween a point and a polygon layer. I can identify wich points intersects my polygon layer but i can't use them in a new layer. I think i have to store the result in a new point layer with a list, but it not seems to be the good argument.

How can i "populate" a point layer with a list like ?

[<QgsGeometry: Point (505555.56033052341081202 6341490.30133426282554865)>, <QgsGeometry: Point (505555.56033052341081202 6341490.30133426282554865)>, <QgsGeometry: Point (505555.56033052341081202 6341490.30133426282554865)>, <QgsGeometry: Point (505555.56033052341081202 6341490.30133426282554865)>]

Here is my code

# Points
myPoints = "C:\pyqgis\mypoints.shp"
layerPoints = QgsVectorLayer(myPoints, '', 'ogr')

# Polygon
myPolygon = "C:\pyqgis\mypolygon.shp"
layerPolygon = QgsVectorLayer(myPolygon, '', 'ogr')

# Intersection
featureList = []
for i in layerPoints.getFeatures():
        for p in layerPolygon.getFeatures():
            if i.geometry().intersects(p.geometry()):
                featureList.append(geom.asWkt())                

# Create layer with result
intersection = QgsFeature()
intersection.setGeometry(QgsGeometry.fromPointXY(featureList))

# Add it to the map
result = iface.addVectorLayer(intersection, "result", "ogr")
1

You can use a memory layer as in following code:

registry = QgsProject.instance()

# Points
myPoints = "/home/zeito/pyqgis_data/mypoints.shp"
layerPoints = QgsVectorLayer(myPoints, 'points', 'ogr')
registry.addMapLayer(layerPoints)

# Polygon
myPolygon = "/home/zeito/pyqgis_data/mypolygon.shp"
layerPolygon = QgsVectorLayer(myPolygon, 'polygon', 'ogr')
registry.addMapLayer(layerPolygon)

# Intersection
featureList = []
for i in layerPoints.getFeatures():
    for p in layerPolygon.getFeatures():
        if i.geometry().intersects(p.geometry()):
            featureList.append(i.geometry())                

epsg = layerPoints.crs().postgisSrid()

# Create layer with result
uri = "Point?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

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

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(featureList)) ]

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

prov.addFeatures(feats)

registry.addMapLayer(mem_layer)

After running it in Python Console of QGIS 3, I got (with my path's layers) following result:

enter image description here

Intersection points are visualized in red; as expected.

  • Great thank you. – GeoGyro Nov 10 at 10:25

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.