Determining the centre of a larger number of arc segments using QGIS

Question

I am using QGIS.

I have a larger data set of circles (different radii) that partially overlap and are intersected with each other (see screenshot below).

I am now looking for a way to determine the centre points of as many circles as possible. So far, this has only been possible approximately, e.g. via different inner buffers.

Is there a way to determine the centre point of a circle using only a partial section of the circle?

It might also be conceivable to simplify the arcs of the circle (line segments) and then form the central perpendicular of these segments - the point of intersection should be very close to the centre of the circle.

Answer 1

First you will need to define some custom functions. Code for the functions, bellow.

from math import acos, dist
from qgis.core import QgsGeometry

@qgsfunction(args='auto', group='Geometry', referenced_columns=[])
def make_line2(points, _, parent):
    points = [point.asPoint() for point in points]
    return QgsGeometry.fromPolylineXY(points)

@qgsfunction(args='auto', group='Geometry', referenced_columns=[])
def angle3(p1, p2, p3, _, parent):
    p1 = p1.asPoint()
    p2 = p2.asPoint()
    p3 = p3.asPoint()
    
    a = dist(p1, p3)
    b = dist(p2, p3)
    c = dist(p1, p2)
    d = acos((b ** 2 + c ** 2 - a ** 2) / (2 * b * c))
    return d

@qgsfunction(args=-1, group='Arrays', referenced_columns=[])
def array_zip(arrays, feature, parent):
    return [list(i) for i in zip(*arrays)]
    
@qgsfunction(args='auto', group='Arrays', referenced_columns=[])
def array_zip2(arrays, feature, parent):
    return [list(i) for i in zip(*arrays)]

@qgsfunction(args='auto', group='Geometry', referenced_columns=[])
def find_center_3p(points, _, parent):
    try:
        p1, p2, p3 = (point.asPoint() for point in points)
    except AttributeError:
        return None
    else:
        x_1 = p1.x()
        y_1 = p1.y()

        x_2 = p2.x()
        y_2 = p2.y()

        x_3 = p3.x()
        y_3 = p3.y()
        
        eq_1 = ((x_1 ** 2 - x_2 ** 2) + (y_1 ** 2 - y_2 ** 2)) / 2
        eq_2 = ((x_2 ** 2 - x_3 ** 2) + (y_2 ** 2 - y_3 ** 2)) / 2
        
        a = (eq_1 - eq_2 * (y_1 - y_2) / (y_2 - y_3)) / ((x_1 - x_2) - (x_2 - x_3) / (y_2 - y_3) * (y_1 - y_2))
        b = (eq_1 - a * (x_1 -x_2)) / (y_1 - y_2)
        
        return QgsGeometry.fromPointXY(QgsPointXY(a, b))

Note: highly recommended to run all this stuff in projected layers. find_center_3p and angle3 use cartesian measurements.

Open the Geometry by expression tool, set the Input layer parameter to your polygon layer, set the Output geometry type to Line. Then, open the expression builder with the button, go to the Function Editor tab and paste this code in a new file (create file with the button) and finally press the Save and Load Functions button. Go back to the Expression tab and paste this expression.

with_variable(
    name:='output',
    value:=with_variable(
        name:='nodes',
        value:=with_variable(
            name:='points',
            value:=geometries_to_array(
                nodes_to_points($geometry)
            ),
            expression:=array_zip2(
                array_filter(
                    array_foreach(
                        array:=generate_series(
                            1,
                            num_points(
                                $geometry
                            )
                        ),
                        expression:=array(
                            @element,
                            180 - degrees(
                                angle3(
                                    @points[@element - 1],
                                    @points[@element],
                                    @points[@element + 1]
                                )
                            )
                        )
                    ),
                    @element[1] > 35 -- threshold for what can be consider a circle intersection
                )
            )[0]
        ),
        expression:=collect_geometries(
            with_variable(
                name:='lines',
                value:=with_variable(
                    name:='vertices',
                    value:=geometries_to_array(
                        nodes_to_points(
                            $geometry
                        )
                    ),
                    expression:=array_foreach(
                        array:=array_cat(
                            array_zip(
                                array_slice(@nodes, 0, -2),
                                array_slice(@nodes, 1, -1)
                            ),
                            array(
                                array(@nodes[-1], -1),
                                array(0, @nodes[0])
                            )
                        ),
                        expression:=make_line2(
                            array_slice(
                                @vertices,
                                @element[0],
                                @element[1]
                            )
                        )
                    )
                ),
                expression:=array_append(
                    array_slice(@lines, 0, -3),
                    combine(@lines[-2], @lines[-1])
                )
            )
        )
    ),
    expression:=if(@output IS NULL, boundary($geometry), @output)
)

what this expression does is to identify where are vertices with very pronounced angles and then, divide the perimeter of the original polygon at this vertices creating line segments which are the circle fragments (arcs). You can set in the expression the minimum angle to be considered as a circle intersection, (default 35) you can decrease it in order to catch more intersections

Convert the result of the above expression to single part geometries with the Multipart to singleparts tool.

Now you have to place three points in this arcs as mentioned by @Vince and solve the circle equation. So open again the Geometry by expression tool and select as the Input layer parameter the single part geometry, the Output geometry type now will be Point, and the expression.

find_center_3p(
    array_foreach(
        generate_series(0, 2),
        line_interpolate_point(@geometry, @element * length($geometry) / 3)
    )
)

Answer 2

1. Points along geometry to create points along the buffer borders
2. Create perpendicular lines for each point with geometry by expression: extend(make_line($geometry, project ($geometry, 1500, radians("angle"-90))), 1500, 0) I use 1500 m to make each line reach buffer centers.
3. Intersect the lines with themselves to produce many points. Most of the points/intersections will be in buffer centers.
4. Cluster them to find the point concentrations. I use min cluster size =10 and max distance = 10
5. Extract the points with "CLUSTER_ID" is not null
6. Dissolve by CLUSTER_ID field to create multipoints at each buffer center
7. Calculate centroid of each multipoint cluster

You will have to play around with points along geometry distance and the cluster parameters for it to work.

---