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")
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)
|improve this answer|||||
  • 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 '18 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 '18 at 9:17
  • @Neven - Honestly, for me the best place to learn was here (GIS:Stack Exchange). – Joseph Nov 28 '18 at 10:58
  • 1
    @Neven - Most welcome and good luck! :) – Joseph Nov 29 '18 at 15:20

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.