5

I have a polygon layer where there may be polygons and circular arcs as in the image.

How can I identify the polygons which have those arcs parts?

enter image description here

8
  • what does geom.wkbType() return?
    – Elio Diaz
    Commented Dec 21, 2021 at 15:52
  • It would be really really hard to identify polygons with such "circle-ish" parts computationally. What is your goal? Maybe a different strategy would be possible. Commented Dec 21, 2021 at 15:55
  • @ElioDiaz it return 6 -Multipolygon
    – Nil
    Commented Dec 21, 2021 at 16:01
  • 4
    Every polygon has segments which could be circular parts, even if it was only two points of a truly massive circle. You need to define your problem in terms of N consecutive vertices which share a common center point when cast as an arc. This is extemely computationally expensive due to all the trig functions, so you can expect this evaluation to slam the CPU on which it is run.
    – Vince
    Commented Dec 21, 2021 at 16:49
  • 4
    I'd approach this rather naive and try to determine a threshold of vertices per area (of bounding box) above which the likeliness of an arc within the polygon is high enough.
    – Erik
    Commented Dec 21, 2021 at 17:03

1 Answer 1

6

Working on the assumption that in an arc like in the posted screenshot, the distance between vertices will be very similar and the azimuth from one vertex to the next will change at a fairly constant rate, I cobbled together this messy code in the Python Console to test the idea:

lyr = iface.activeLayer()

for f in lyr.selectedFeatures():
    geom = f.geometry()

    for p in geom.parts():
        verts = p.vertices()
        p0 = verts.next()
        p1 = verts.next()
        dist = p0.distance(p1)
        az = p0.azimuth(p1)
        dev = 0
        tol = 0.1
        angular_tol_max = 10
        angular_tol_min = 0.01
        in_arc = False
        
        for v in verts:
            d = p1.distance(v)
            a = p1.azimuth(v)
            if((abs(d-dist)/dist) > tol):
                dist = d
                in_arc = False
            else:
                if (angular_tol_min < abs(az-a) < angular_tol_max) and dev != 0:
                    #print ("this might be an arc",((az-a)/dev))
                    if( (-tol) < ((az-a)/dev) > (tol)):     
                        #print ("this realy might be an arc")
                        if in_arc == False:
                            rb = QgsRubberBand(iface.mapCanvas())
                            rb.setColor(QColor(255, 0, 0))
                            rb.setWidth(3)
                            rb.addPoint(QgsPointXY(p0.x(), p0.y()))
                            rb.addPoint(QgsPointXY(p1.x(), p1.y()))
                        rb.addPoint(QgsPointXY(v.x(), v.y()))
                        in_arc = True
                    else:
                        in_arc = False
                else:
                    in_arc = False
            p0 = p1
            p1 = v
            dev = az - a
            az = a

It iterates through the vertices comparing the distance and azimuth to the previous vertex, when it detects a segment with a similar length and change in direction as the previous segment it assumes they are part of the arc. The tolerance can be changed by the value of the tol variable (current value of 0.1 means up to 10% difference in segment length) and the angular_tol_max and angular_tol_min currently 10 and 0.01 meaning the azimuth change must be at less than 10 degrees and at at least 0.01 degrees.

It seems to work OK, I expect it would fail where the azimuth of adjacent segments changes from slightly less than 360 to slightly more than zero.

enter image description here

2
  • It works i.imgur.com/MxzaQ1J.png Thanks a lot!
    – Nil
    Commented Dec 22, 2021 at 7:09
  • 1
    You should publish this algorithm and write it in a fast implementation. It would get used. Cheers,
    – dgketchum
    Commented Nov 10, 2022 at 2:03

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.