PyQGIS Polyline Creation – Creating Perpendicular Line between Two Points Inside Polygon Using PyQGIS

perpendicularpoint-in-polygonpolyline-creationpyqgis

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)

Best Answer

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:

Related Solutions

[GIS] Offsetting points perpendicular to line using ArcMap

I found a way! I used a couple of different sources and manipulated the script to fit my needs (mostly from gotanuki on StackExchange). A variation of this is already on the site, but the of data that I was using is different. In case someone else runs into this issue, this is the solution I used:

# Define orientation (start, mid, end) and draw Matchlines function
# This function is defined before it is called
def orientMatchlines(feat,ix,iy):
        # If the line is horizontal or vertical, the slope and
        # negative reciprocal calculations will fail, so do this instead
        if starty==endy or startx==endx:
            if starty == endy:
                y1 = iy + matchDist
                y2 = iy - matchDist
                x1 = ix
                x2 = ix
            if startx == endx:
                y1 = iy
                y2 = iy
                x1 = ix + matchDist
                x2 = ix - matchDist
        else:
            # Get slope of line
            m = float((starty - endy)/(startx - endx))

            # Get negative reciprocal
            negativereciprocal = -1*((startx - endx)/(starty - endy))

            # For all values of slope, calculate perpendicular line
            # with length = matchDist
            if m > 0:
                if m >= 1:
                    y1 = negativereciprocal*(matchDist)+ iy
                    y2 = negativereciprocal*(-matchDist) + iy
                    x1 = ix + matchDist
                    x2 = ix - matchDist
                if m < 1:
                    y1 = iy + matchDist
                    y2 = iy - matchDist
                    x1 = (matchDist/negativereciprocal) + ix
                    x2 = (-matchDist/negativereciprocal)+ ix
            if m < 0:
                if m >= -1:
                    y1 = iy + matchDist
                    y2 = iy - matchDist
                    x1 = (matchDist/negativereciprocal) + ix
                    x2 = (-matchDist/negativereciprocal)+ ix
                if m < -1:
                    y1 = negativereciprocal*(matchDist)+ iy
                    y2 = negativereciprocal*(-matchDist) + iy
                    x1 = ix + matchDist
                    x2 = ix - matchDist
        point1.X = x1
        point1.Y = y1
        point2.X = x2
        point2.Y = y2
        lineArray.add(point1)
        lineArray.add(point2)

        del x1
        del x2
        del y1
        del y2

# Create search cursor in inLine
rows = arcpy.SearchCursor(inLine)

# Get number of records in inLine
numRecords = int(arcpy.GetCount_management(inLine).getOutput(0))

# Create new point files and array to collect values
point1 = arcpy.Point()
point2 = arcpy.Point()
lineArray = arcpy.Array()

# Define counter
counter = 0

# Loop over rows in outLine
for row in rows:
    # Create the geometry object
    feat = row.Shape

    # Get coordinate values as lists
    firstpoint = feat.firstPoint
    lastpoint = feat.lastPoint
    midpoint = feat.centroid

    # Get X and Y values for each point
    startx = firstpoint.X
    starty = firstpoint.Y
    endx = lastpoint.X
    endy = lastpoint.Y
    midx = midpoint.X
    midy = midpoint.Y

    m = ((starty - endy)/(startx - endx+.0001))

    # For all points besides the last one
    if counter < numRecords - 1:
        orientMatchlines(feat,startx,starty)
    # For the last point
    else:
        orientMatchlines(feat,endx,endy)

    #Create insert cursor
    cur = arcpy.InsertCursor(matchlines)

    #Insert new row from array
    feat = cur.newRow()
    #feat.slope = m
    feat.shape = lineArray
    cur.insertRow(feat)
    lineArray.removeAll()

    del cur

    # Increase counter by 1 and start again
    counter = counter + 1

del row
del rows

[GIS] Draw a perpendicular line between point and line layer

You need a reference to the layers. I do that with QgsMapCanvas class in the next code. PyQGIS also has classes to find 'Closest Segments'. For this reason you can avoid to use your 'intersect_point_to_line' function. I used 'closestSegmentWithContext' of QgsGeometry instead.

mapcanvas = iface.mapCanvas()

layers = mapcanvas.layers()

feat_points = [ feat for feat in layers[0].getFeatures() ]

feat_line = layers[1].getFeatures().next()

geoms = [ feat_line.geometry().closestSegmentWithContext(feat.geometry().asPoint())
          for feat in feat_points ]

#geom is a tupla. I need only second term, geom[1], that is a QgsPoint
for i, geom in enumerate(geoms):
    closest_line = QgsGeometry.fromPolyline([ feat_points[i].geometry().asPoint(), 
                                             geom[1] ])

    print closest_line.exportToWkt()

I employed above code with next layers; where you can have two possible situations.

After running the code, at the Python Console of QGIS were printed two Line Strings in WKT format. By using QuickWKT plugin of QGIS I got:

You can observe that, in one case, second closest segment is not perpendicular.

Best Answer

Related Solutions

[GIS] Offsetting points perpendicular to line using ArcMap

[GIS] Draw a perpendicular line between point and line layer

Related Question