1

I want to extract points which are overlapped by polygons to a shapefile with PyQGIS. I typed the code below, I don't get an error, but I get points that I have on the beginning. So my intersection doesn't work, and I can't see why.

points = QgsVectorLayer(points_path, "layer", "ogr")  
polygons = QgsVectorLayer(polygon_path, "layerpolygon", "ogr") 
int_points = QgsVectorLayer("Point", "inttersected_points", "memory")
pr = int_points.dataProvider()

for polygon in polygons.getFeatures():  
    geom_polygon = polygon.geometry()  
    features_point = points.getFeatures(QgsFeatureRequest().setFilterRect(geom_polygon.boundingBox()))
    for one_feature_point in features_point:

        if one_feature_point.geometry().intersected(geom_polygon):
             int_points = one_feature_point.getFeatures() 

QgsVectorFileWriter.writeAsVectorFormat(int_points, save_int_as_shp, "CP120", points.crs(), "ESRI Shapefile")

1 Answer 1

2

I personally find it easier to find intersected points and store them in a list. Once all points have been found, you could then select them all and save them in a shapefile in one go. So you could try something like the following:

points = QgsVectorLayer(points_path, "layer", "ogr")  
polygons = QgsVectorLayer(polygon_path, "layerpolygon", "ogr") 

# Create list to store id of intersected points
ids = []

for polygon in polygons.getFeatures():
    geom_polygon = polygon.geometry()
    features_point = points.getFeatures(QgsFeatureRequest().setFilterRect(geom_polygon.boundingBox()))
    for one_feature_point in features_point:
        if one_feature_point.geometry().intersects(geom_polygon):
            # Store id of intersected point
            ids.append(one_feature_point.id())

# Select intersected points
points.selectByIds(ids)
# Write to output
QgsVectorFileWriter.writeAsVectorFormat(points, output_path, "CP120", points.crs(), "ESRI Shapefile", onlySelected=True)
4
  • For some reason, I am getting empty shapefiles on the output. But I know that there is a intersect because I do that intersect manually from a few examples.
    – Neven
    Nov 23, 2018 at 9:14
  • Btw, thank you for the reply! Your code is intuitive and great for understanding! Can you please give me the advice about how to learn all the classes and etc for PyQGIS. A read Cookbook, but it's a for 2.0 API.
    – Neven
    Nov 23, 2018 at 9:17
  • @Neven - Honestly, for me the best place to learn was here (GIS:Stack Exchange).
    – Joseph
    Nov 28, 2018 at 10:58
  • 1
    @Neven - Most welcome and good luck! :)
    – Joseph
    Nov 29, 2018 at 15:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.