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.

  • 1
    en.wikipedia.org/wiki/Rotating_calipers is something like this you are looking for? – Elio Diaz Jun 24 at 16:49
  • what's the distance the plane needs to take off? – Elio Diaz Jun 24 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. – Michael Sousa Jun 25 at 3:42
  • The correct answer should become a plugin! – Stu Smith Jun 26 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'.

| improve this answer | |
  • 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 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 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? – Michael Sousa Jun 25 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 at 15:23
4

When holes of the polygons have to be considered

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

input

Step 1. Proceed with "Delete holes".

step_1

Step 2. Apply "Polygons to lines"

step_2

Step 3. Use "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 "Join by lines (hub lines)". Afterwards the application of "Fix geometries", "Remove null geometries" and "Delete duplicate geometries" is probable.

step_4

Step 5. Proceed with "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 "Extract by location" (are within for geometrical predicate) for the Result of the Step 5.

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

step_7

| improve this answer | |
3

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

| improve this answer | |

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.