4

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.

enter image description here

3
  • 2
    Please decide which software you want to use.
    – Erik
    Oct 26, 2022 at 9:03
  • 1
    Welcome to GIS SE! As a new user be sure to take the Tour to learn about our focussed Q&A format. Please decide which of QGIS and ArcGIS Pro you wish to ask about within this question so that it and its answers can be focused on what you have tried and where you are stuck.
    – PolyGeo
    Oct 26, 2022 at 9:04
  • 2
    It is certainly possible to identify the center if you can identify at least three points along the curve of that circle -- just construct a chord between the two outer points, bisect it, then extend the line from the middle point through the midpoint until the lines from all three points are equal in length.
    – Vince
    Oct 26, 2022 at 12:37

2 Answers 2

3

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 2 button, go to the Function Editor tab and paste this code in a new file (create file with the 4 button) and finally press the Save and Load Functions 5 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)
    )
)

7

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.

enter image description here


enter image description here

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.