3

I have a polygon shp with agricultural fields and a point shp representing GPS locations. I would like to create separate point layers for each field the features intersect. Below is a screenshot of my data.

enter image description here

A selection the code would need to make and then extract looks roughly like this:

enter image description here

I thought of iterating through each polygon feature and extracting all points that intersect with it. This question aims for exactly that, but the solution also extracts empty intersections, giving me an insane amount of useless outputs which also freeze QGIS. Here is the code given in the solution:

input_shp = QgsProject.instance().mapLayersByName('shp_1')[0]
intersect_shp = QgsProject.instance().mapLayersByName('shp_2')[0]

for feature in intersect_shp.getFeatures():
    intersect_shp.selectByIds([feature.id()])
    processing.runAndLoadResults("native:extractbylocation",
    {'INPUT':input_shp,
    'PREDICATE':[6],
    'INTERSECT':QgsProcessingFeatureSourceDefinition(intersect_shp.id(),
    selectedFeaturesOnly=True,
    featureLimit=-1,
    geometryCheck=QgsFeatureRequest.GeometryAbortOnInvalid),
    'OUTPUT':'TEMPORARY_OUTPUT'})

intersect_shp.selectByIds([])

I tried adding an "if" statement to skip empty features, but they are still being loaded.

input_shp = QgsProject.instance().mapLayersByName('shp_1')[0]
intersect_shp = QgsProject.instance().mapLayersByName('shp_2')[0]

for feature in intersect_shp.getFeatures():
    if not feature.geometry():
        continue
    intersect_shp.selectByIds([feature.id()])
    processing.runAndLoadResults("native:extractbylocation",
    {'INPUT':input_shp,
    'PREDICATE':[6],
    'INTERSECT':QgsProcessingFeatureSourceDefinition(intersect_shp.id(),
    selectedFeaturesOnly=True,
    featureLimit=-1,
    geometryCheck=QgsFeatureRequest.GeometryAbortOnInvalid),
    'OUTPUT':'TEMPORARY_OUTPUT'})

intersect_shp.selectByIds([])

My workaround at the moment is to let the original code run normally, but delete the empty layers in the end before adding them to the map. However, this is still intensive for QGIS (it still freezes for a moment) and probably not good for bigger datasets.

input_shp = QgsProject.instance().mapLayersByName('shp_1')[0]
intersect_shp = QgsProject.instance().mapLayersByName('shp_2')[0]

for feature in intersect_shp.getFeatures():
    intersect_shp.selectByIds([feature.id()])
    result = processing.runAndLoadResults("native:extractbylocation",
    {'INPUT':input_shp,
    'PREDICATE':[6],
    'INTERSECT':QgsProcessingFeatureSourceDefinition(intersect_shp.id(),
    selectedFeaturesOnly=True,
    featureLimit=-1,
    geometryCheck=QgsFeatureRequest.GeometryAbortOnInvalid),
    'OUTPUT':'TEMPORARY_OUTPUT'})
    output_layer_name = result['OUTPUT']
    output_layer = QgsProject.instance().mapLayer(output_layer_name)
    if not output_layer.featureCount():
        QgsProject.instance().removeMapLayer(output_layer.id())
    else:
        QgsProject.instance().addMapLayer(output_layer)

intersect_shp.selectByIds([])

As such, I'd like to use the for-loop to avoid empty features altogether. But why is that code not working?

I have checked my polygon layer for validity and null geometries. There are no problems.

2
  • 1
    I'm running your last chunk of code with the points layer as input_shp and polygon layer as intersects_shp and I am not getting any 'empty' layers. It is behaving as expected.
    – Cushen
    Commented May 12, 2023 at 7:47
  • If the last chunk is the workaround, it does work, it's just not ideal. It is the second mentioned code that still leaves empty layers up. Commented May 15, 2023 at 14:33

1 Answer 1

4

I suggest you run join attributes by location first, in order to add the ID or your parcels to your points, which allows you to group/categorize your point data.

Then run split vector layer on the joined points, using the joined parcel-ID as the unique key to split by.

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.