8

I'm looking for an algorithm to draw max distance straight lines contained in a polygon (without crossing polygon boundaries). How can I perform it?

As you can see some polygons are very oddly shaped.

enter image description here

The area I am working with is millions of acres and contains hundreds of thousands of individual polygons. From what I can tell the polygons are pretty clean (no holes, overlaps, etc.)

The practical use of this is to filter what bodies of water are capable of having a floatplane land and take off on them.

The only criteria I have for a solution is that I can perform it in QGIS and that the lines are straight.

4
  • 1
    en.wikipedia.org/wiki/Rotating_calipers is something like this you are looking for?
    – Elio Diaz
    Jun 24 '20 at 16:49
  • what's the distance the plane needs to take off?
    – Elio Diaz
    Jun 24 '20 at 21:03
  • 1
    I'm not sure how to share the data here but I can email anyone the shapefile if necessary or if they thing it will be helpful. The inner lines can be connected from any point to any point, they don't necessarily need to be vertex to vertex. The distance for the plane is ~1500yds. Jun 25 '20 at 3:42
  • The correct answer should become a plugin!
    – Stu Smith
    Jun 26 '20 at 4:41
7

When holes of the polygons have to be avoided

So, this is an extension to my previous answer Calculating the longest distance within polygon in QGIS but with some changes in the Step 3, particularly in the query.

SELECT p1.id, setsrid(make_line(p1.geometry, p2.geometry),  #put your srid here),
       max(st_length(make_line(p1.geometry, p2.geometry))) AS length
FROM "Points" AS p1, "polygons" AS p
JOIN "Points" AS p2 ON p1.id = p2.id
WHERE NOT st_equals(p1.geometry, p2.geometry)
      AND st_within(make_line(p1.geometry, p2.geometry), st_buffer(p.geometry, 0.00005))
GROUP BY p1.id

Note that in the query above additionally the geometry of the original polygons were used.

To be more example-realisting I considered different polygons to those that I had in my previous answer, see image below.

input

The corresponding result will be looking as

result

Note that the result is approximate because a bigger distance was used on the step 'Points along geometry'.

I made with 'Points along geometry', however it can also be done with the result of 'Extract vertices'.

4
  • 1
    This approach is not quite correct for non convex polygons. In general, the lines joining the vertices need to be extended outwards towards the overall polygon bounding box before clipping them with the original polygon. Otherwise you'll get misleading results for shapes such as a "C" or polygons with holes.
    – ndawson
    Jun 25 '20 at 18:31
  • 2
    Actually your demo image shows this situation -- the rightmost polygon can clearly have a longer line drawn within it
    – ndawson
    Jun 25 '20 at 18:33
  • I can get the virtual layer created I just can't get the lines to visualize. They exist in the attribute table and have distances that seem reasonable. Any suggestions? Jun 25 '20 at 20:09
  • You also need to test every vertex in all rings against every other vertex in all rings, extend these lines out, clip to the polygon and then take the longest part returned from the clip result..
    – ndawson
    Jun 26 '20 at 15:23
5

This is also an approximation. It finds the longest line between existing polygon vertices within each polygon. So the more vertices you have the better the results should be (but the executing time will be longer). For example add more vertices using Densify by interval if the results are not good enough. But from what I can see you seem to have alot of vertices.

It can be slow, so try it on a subset of polygons first.

import itertools

layer = iface.activeLayer() #Click layer in tree

#Create empty line layer
vl = QgsVectorLayer("LineString?crs={}&index=yes".format(layer.crs().authid()), "Longest_line", "memory")
provider = vl.dataProvider()

#For each polygon find the longest line that is within the polygon
for feat in layer.getFeatures():
    verts = [v for v in feat.geometry().vertices()] #List all vertices
    all_lines = []
    for p1,p2 in itertools.combinations(verts, 2): #For every combination of two vertices
        all_lines.append(QgsGeometry.fromPolyline([p1,p2])) #Create a line
    all_lines = [line for line in all_lines if line.within(feat.geometry())] #Check if line is within polygon
    if len(all_lines)>0:
        longest_line = max(all_lines, key=lambda x: x.length()) #Find longest line
        #Create a line feature from the longest line within polygon
        f = QgsFeature()
        f.setGeometry(longest_line)
        provider.addFeature(f)

QgsProject.instance().addMapLayer(vl)

enter image description here

0
4

When holes of the polygons have to be considered

Let's assume there is a polygon layer 'polygons' with its corresponding attribute table accordingly, see image below.

input

Step 1. Proceed with the "Delete holes"

step_1

Step 2. Apply the "Polygons to lines"

step_2

Step 3. Use the "Points along geometry". The output of this algorithm save additionally as a permanent file. Both layers will be used at the Step 4

step_3

Step 4. Make use of the "Join by lines (hub lines)". Afterwards the application of the "Fix geometries", "Remove null geometries" and "Delete duplicate geometries" is probable

step_4

Step 5. Proceed with the "Clip" between the result of the Step 4 and initial polygons

step_5

Step 6. Apply a tiny "Buffer" for initial polygons. And after make use of the "Extract by location" (are within for geometrical predicate) for the Result of the Step 5

Step 7. Use the "Extract by expression" using the following expression $length = maximum($length, "id").

step_7

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.