Creating perpendicular line between two points inside polygon using PyQGIS

Question

I'm looking for a solution in PyQGIS to create a cut line at the midpoint of another line that connects two points within a polygon. The objective is to cut the polygon into two parts.

In fact, this would be the alternative solution to the script below used to generate subdivisions in polygons with more than 2 points on it (based on cluster points and Voronoi). My problem is when I have two points inside one of these polygons.

The code was adapted from this link.

def split_polygon(poligonfeatures, numberofpoints, numberofparts):
    # Random points inside polygons
    alg_process = 'RandomPointsInsidePolygons'
    alg_params = {
        'INPUT': poligonfeatures,
        'MIN_DISTANCE': None,
        'STRATEGY': 0,
        'VALUE': numberofpoints,
        'OUTPUT': 'TEMPORARY_OUTPUT'
    }
    outputs[alg_process] = processing.run('qgis:randompointsinsidepolygons', alg_params)

    # K-means clustering
    alg_process = 'KmeansClustering'
    alg_params = {
        'CLUSTERS': numberofparts,
        'FIELD_NAME': 'CLUSTER_ID',
        'INPUT': outputs['RandomPointsInsidePolygons']['OUTPUT'],
        'OUTPUT': 'TEMPORARY_OUTPUT'
    }
    outputs[alg_process] = processing.run('native:kmeansclustering', alg_params)

    # Concave hull (k-nearest neighbor)
    alg_process = 'ConcaveHullKnearestNeighbor'
    alg_params = {
        'FIELD': 'CLUSTER_ID',
        'INPUT': outputs['KmeansClustering']['OUTPUT'],
        'KNEIGHBORS': 3,
        'OUTPUT': 'TEMPORARY_OUTPUT'
    }
    outputs[alg_process] = processing.run('qgis:knearestconcavehull', alg_params)

    # Centroids
    alg_process = 'Centroids'
    alg_params = {
        'ALL_PARTS': False,
        'INPUT': outputs['ConcaveHullKnearestNeighbor']['OUTPUT'],
        'OUTPUT': 'TEMPORARY_OUTPUT'
    }
    outputs[alg_process] = processing.run('native:centroids', alg_params)

    # Voronoi polygons
    alg_process = 'VoronoiPolygons'
    alg_params = {
        'BUFFER': 100,
        'INPUT': outputs['Centroids']['OUTPUT'],
        'OUTPUT': 'TEMPORARY_OUTPUT'
    }
    outputs[alg_process] = processing.run('qgis:voronoipolygons', alg_params)

    # Intersection
    alg_process = 'Intersection'
    alg_params = {
        'INPUT': poligonfeatures,
        'INPUT_FIELDS': [''],
        'OVERLAY': outputs['VoronoiPolygons']['OUTPUT'],
        'OVERLAY_FIELDS': [''],
        'OVERLAY_FIELDS_PREFIX': '',
        'OUTPUT': 'TEMPORARY_OUTPUT'
    }
    outputs[alg_process] = processing.runAndLoadResults('native:intersection', alg_params)

split_polygon(layer, 50, 2)

Answer 1

In this solution, none of the algorithms from the processing was used.

Let's assume there is a polygon layer called 'polygon' (five features) and a point layer called 'points_in_polygon' (12 features) with its attribute tables, see the image below.

Proceed with Plugins > Python Console > Show Editor and paste the script below:

# imports
from PyQt5.QtCore import QVariant
from qgis.core import QgsProject, QgsField, QgsFeature, QgsVectorLayer

# layers' names
point_layer = QgsProject.instance().mapLayersByName('points_in_polygon')[0]
polygon_layer = QgsProject.instance().mapLayersByName('polygon')[0]

# crs for the output
output_crs = polygon_layer.crs().authid()

# creating a polyline layer for the output
polyline_layer = QgsVectorLayer(f"LineString?crs={output_crs}", "temp", "memory")
provider = polyline_layer.dataProvider()
provider.addAttributes([QgsField("id", QVariant.Int)])
polyline_layer.updateFields()

# looping over polygon and point features
polygons_with_points = {}

for polygon in polygon_layer.getFeatures():
    polygon_geom = polygon.geometry()
    points_in_polygon = []
    
    for point in point_layer.getFeatures():
        point_geom = point.geometry()
        if point_geom.within(polygon_geom):
            points_in_polygon.append(point_geom.asPoint())

    # only polygons with exactly two points inside are needed
    if len(points_in_polygon) == 2:
        polygons_with_points[polygon.id()] = QgsGeometry.fromPolylineXY(points_in_polygon)

# creating and editing polyline features
polyline_layer.startEditing()
for id, geom in polygons_with_points.items():
    feat = QgsFeature(polyline_layer.fields())
    feat.setAttribute("id", id)
    # getting line centroid as a rotation center
    centroid = geom.centroid().asPoint()
    # rotating line clockwise in 90 degrees
    geom.rotate(90, centroid)
    # extending line from start and end in 20 meters
    new_geom = geom.extendLine(0.0, float(20)).extendLine(float(20), 0.0)
    feat.setGeometry(new_geom)
    polyline_layer.addFeature(feat)
polyline_layer.endEditCommand()
polyline_layer.commitChanges()

# adding output layer to the canvas
QgsProject.instance().addMapLayer(polyline_layer)

Change the names of the input layer in the first lines of code. Press Run script and get the output that will look like this:

Note: For line extension, probably a more robust technique should be used.

---

References:

• Creating a Line feature in PyQGIS QGIS 3.x
• How can I define the layer CRS for a memory layer in pyqgis?
• Creating polygon from list of coordinates: QgsGeometry.fromPolygon has unexpected type list
• What's the difference between QgsPoint, QgsPointXY and QgsGeometry.fromPointXY(QgsPointXY()) in PyQGIS?
• Shortening line with PyQGIS
• Extending only one end of line using PyQGIS

Answer 2

If it is still necessary to use the code adapted from QGISPolygonSplitting.py, there are several suggestions to make it work:

• The holder for output is missing i.e. outputs = {} at the beginning

• parameter numberofparts must be greater than 2 i.e. numberofparts > 2, otherwise the QgsProcessingException will raise when applying the "qgis:voronoipolygons":

  results, ok = alg.run(parameters, context, feedback, {}, False)
  _core.QgsProcessingException: Input file should contain at least 3 points.
  Choose another file and try again.
  

So, after all your code will look as follows:

def split_polygon(poligonfeatures, numberofpoints, numberofparts):

    outputs = {}

    # Random points inside polygons
    alg_process = 'RandomPointsInsidePolygons'
    alg_params = {
        'INPUT': poligonfeatures,
        'MIN_DISTANCE': None,
        'STRATEGY': 0,
        'VALUE': numberofpoints,
        'OUTPUT': 'TEMPORARY_OUTPUT'
    }
    outputs[alg_process] = processing.run('qgis:randompointsinsidepolygons', alg_params)

    # K-means clustering
    alg_process = 'KmeansClustering'
    alg_params = {
        'CLUSTERS': numberofparts,
        'FIELD_NAME': 'CLUSTER_ID',
        'INPUT': outputs['RandomPointsInsidePolygons']['OUTPUT'],
        'OUTPUT': 'TEMPORARY_OUTPUT'
    }
    outputs[alg_process] = processing.run('native:kmeansclustering', alg_params)

    # Concave hull (k-nearest neighbor)
    alg_process = 'ConcaveHullKnearestNeighbor'
    alg_params = {
        'FIELD': 'CLUSTER_ID',
        'INPUT': outputs['KmeansClustering']['OUTPUT'],
        'KNEIGHBORS': 3,
        'OUTPUT': 'TEMPORARY_OUTPUT'
    }
    outputs[alg_process] = processing.run('qgis:knearestconcavehull', alg_params)

    # Centroids
    alg_process = 'Centroids'
    alg_params = {
        'ALL_PARTS': False,
        'INPUT': outputs['ConcaveHullKnearestNeighbor']['OUTPUT'],
        'OUTPUT': 'TEMPORARY_OUTPUT'
    }
    outputs[alg_process] = processing.run('native:centroids', alg_params)

    # Voronoi polygons
    alg_process = 'VoronoiPolygons'
    alg_params = {
        'BUFFER': 100,
        'INPUT': outputs['Centroids']['OUTPUT'],
        'OUTPUT': 'TEMPORARY_OUTPUT'
    }
    outputs[alg_process] = processing.run('qgis:voronoipolygons', alg_params)

    # Intersection
    alg_process = 'Intersection'
    alg_params = {
        'INPUT': poligonfeatures,
        'INPUT_FIELDS': [''],
        'OVERLAY': outputs['VoronoiPolygons']['OUTPUT'],
        'OVERLAY_FIELDS': [''],
        'OVERLAY_FIELDS_PREFIX': '',
        'OUTPUT': 'TEMPORARY_OUTPUT'
    }
    outputs[alg_process] = processing.runAndLoadResults('native:intersection', alg_params)

    
split_polygon('polygon', 50, 3)