Removing points outside of polygon but keeping points with the same feature ID using PyQGIS

Question

I have a polygon layer and a point layer. The point layer has specific features in its table such as an ID called "fwg_type". There can be several points with the same "fw_type".

I want to remove all the points that lie outside the polygon except those that share the same ID as points, which lie inside the layer boundaries.

So if points with ID 1 lie inside and outside the layer, they would not be removed.

Here is my script. It does not work yet. Any ideas?

# Replace 'point_layer' and 'polygon_layer' with the actual names of your layers
point_layer = QgsProject.instance().mapLayersByName('point_layer')[0]
polygon_layer = QgsProject.instance().mapLayersByName('polygon_layer')[0]

# Create a set to store the IDs of points within the polygon
ids_within_polygon = set()

# Iterate through the features in the polygon layer and store their IDs
for feature in polygon_layer.getFeatures():
    ids_within_polygon.add(feature['ID'])  # Replace 'ID' with the actual field name

# Create a new memory layer to store the filtered features
new_layer = QgsVectorLayer("Point?crs=" + point_layer.crs().authid(), "Filtered Points", "memory")
new_layer.startEditing()

# Define the attributes and fields for the new layer
new_layer.dataProvider().addAttributes(point_layer.fields())
new_layer.updateFields()

# Create a list to store the features to keep
features_to_keep = []

# Iterate through the features in the point layer
for feature in point_layer.getFeatures():
    if feature['ID'] in ids_within_polygon or feature['ID'] not in ids_within_polygon:
        features_to_keep.append(feature)

# Add the features to keep to the new layer
new_layer.addFeatures(features_to_keep)
new_layer.commitChanges()

# Add the new layer to the QGIS map
QgsProject.instance().addMapLayer(new_layer)

Answer 1

There are some logic omissions and errors in the posted script. In addition to not including any spatial geometric testing between points and polygons, the ids_within_polygon variable is storing "ID" values of polygon features not points and this conditional test in the point feature iteration

if feature['ID'] in ids_within_polygon or feature['ID'] not in ids_within_polygon:

contains two opposite statements separated by an or operator; therefore every feature will satisfy one or the other statement.

Try the following. You will need check the layer names for your point and polygon layers as well as the name of the relevant ID field in your point layer (in your question you mentioned both "fwg_type" and "fw_type"). In this answer I used "fwg_type" so make sure this is correct.

from qgis.core import QgsProject, QgsVectorLayer, QgsGeometry, QgsCoordinateTransform

project = QgsProject.instance()

point_lyr = project.mapLayersByName('points')[0]

polygon_lyr = project.mapLayersByName('polygons')[0]

output_lyr = QgsVectorLayer(f'Point?crs={point_lyr.crs().authid()}', 'result', 'memory')

output_lyr.dataProvider().addAttributes(point_lyr.fields())

output_lyr.updateFields()

poly_combined_geom = QgsGeometry.collectGeometry([ft.geometry() for ft in polygon_lyr.getFeatures()])

# transform combined polygon geometry to point layer crs if required
if point_lyr.crs() != polygon_lyr.crs():
    xform = QgsCoordinateTransform(polygon_lyr.crs(), point_lyr.crs(), project)
    poly_combined_geom.transform(xform)

geom_engine = QgsGeometry.createGeometryEngine(poly_combined_geom.constGet())

geom_engine.prepareGeometry()

pts_in = [ft for ft in point_lyr.getFeatures() if geom_engine.contains(ft.geometry().constGet())]

ids_in = set(ft['fwg_type'] for ft in pts_in)

pts_out_to_keep = [ft for ft in point_lyr.getFeatures() if geom_engine.disjoint(ft.geometry().constGet()) and ft['fwg_type'] in ids_in]

all_fts = pts_in + pts_out_to_keep

output_lyr.dataProvider().addFeatures(all_fts)

output_lyr.updateExtents()

project.addMapLayer(output_lyr)

Input layers (points labelled with "fwg_type" field):

Result:

---

References:

• QGIS API Documentation | QgsGeometryEngine class

Answer 2

It works

I tried adding another point layer and keep its points if the first point layer has points within the polygon with the same value.

I added the following to the code but it does not seem to work. At least he just prints the remaining points of the first point layer.

from qgis.core import QgsProject, QgsVectorLayer, QgsGeometry, QgsCoordinateTransform

project = QgsProject.instance()

point_lyr = project.mapLayersByName('05334036_STR_FPG')[0]
point_lyr_2 = project.mapLayersByName('05334036_STR_WGF')[0]
polygon_lyr = project.mapLayersByName('L-den / dB(A)')[0]

output_lyr = QgsVectorLayer(f'Point?crs={point_lyr.crs().authid()}', 'result', 'memory')

output_lyr.dataProvider().addAttributes(point_lyr.fields())
output_lyr.dataProvider().addAttributes(point_lyr_2.fields())
output_lyr.updateFields()

poly_combined_geom = QgsGeometry.collectGeometry([ft.geometry() for ft in polygon_lyr.getFeatures()])

# transform combined polygon geometry to point layer crs if required
if point_lyr.crs() != polygon_lyr.crs():
    xform = QgsCoordinateTransform(polygon_lyr.crs(), point_lyr.crs(), project)
    poly_combined_geom.transform(xform)

geom_engine = QgsGeometry.createGeometryEngine(poly_combined_geom.constGet())

geom_engine.prepareGeometry()

pts_in = [ft for ft in point_lyr.getFeatures() if geom_engine.contains(ft.geometry().constGet())]
pts_in_2 = [ft for ft in point_lyr_2.getFeatures() if geom_engine.contains(ft.geometry().constGet())]


ids_in = set(ft['GEB_ID'] for ft in pts_in)
ids_in_2 = set(ft['GEB_ID'] for ft in pts_in_2)

pts_out_to_keep = [ft for ft in point_lyr.getFeatures() if geom_engine.disjoint(ft.geometry().constGet()) and ft['GEB_ID'] in ids_in]
pts_out_to_keep_2 = [ft for ft in point_lyr_2.getFeatures() if geom_engine.disjoint(ft.geometry().constGet()) and ft['GEB_ID'] in ids_in]

all_fts = pts_in + pts_in + pts_out_to_keep + pts_out_to_keep_2

output_lyr.dataProvider().addFeatures(all_fts)

output_lyr.updateExtents()

project.addMapLayer(output_lyr)