4

I'm working with QGIS 3.16 Hannover. I've got a point layer (shapefile) with 3,000 features and I want to create the smallest convex hull that contains 90% of those points. Well, ideally, I would like to obtain an onion peeling set of minimum-area polygons, containing an ascending percentage of points of the layer.

Is it possible to do so?

2
  • Like an alpha-trimmed mean convex polygon?
    – wingnut
    Apr 25 at 14:09
  • 7
    I don't think it's possible to do it deterministically in reasonable time. Too many combinations, and it appears np-complete. Sounds a bit like an alpha trimmed mean. Maybe start with a convex polygon for all points, then delete the point in the polygon farthest from the centroid, and iterate. No guarantees.
    – wingnut
    Apr 25 at 14:18
8

The basic idea The basic idea of this solution is to find the "most central" point of all your points - meaning the point with the shortest accumulated distance to all other (or the 90% closest) points. The idea behind this: the shorter the distance, the smaller the area needed to cover all of them. From this central point than select the 90% (2700) that are closest to it and create a convex hull to these points.

This method is a quick way to get close to a solution. Probably the only way to get a 100% accurate solution is to create convex hulls for all 3000 points separately, calculating the area and keeping the smallest one. This uses extreme ressources and is probably not realistic. See some thoughts abbout accuracy of this solution at the bottom.

For the "onion peeling set of minimun-area polygons" see below - it's easy to create once you have implemented the basic solution.

In more details:

Preparation I created 3000 random points with an extent of the layer of ca. 170 x 130 km. This layer is called 3000_points (change it accordingly in step 3 if you use another name).

Implementation:

  1. I created a new attribute that calculates for every point the distance to the 2700 (90 % of 3000) closest points and than sums this up - thus aggregates the values (pseudocode): distance to closest + distance to second closest + distance to third closest + (...) + distance to 2700th closest point. I used this expression, it took around three minutes for caluculation. Be aware that overlay_nearest() function is available since QGIS 3.16, array_sum () is available since QGIS 3.18.
array_sum ( 
    array_foreach (
        overlay_nearest( 
            @layer, 
            $geometry, 
            limit:=2700
        ),
        length (
            make_line (
                $geometry, 
                @element
            )
        )
    )
)
  1. Now sort the attribute table to find the central point: it's the one with the smallest value in the attribute created. This point automatically will be close to the center of the layer as only from there, distance to (most) other points tends to be low. Copy this point to a new layer.

  2. Now from this new layer, containing only one feature (the center point), again use overlay_nearest( ) to find the 2700 closest points on the orginal layer 3000_points:

overlay_nearest (
    '3000_points',
    $geometry,
    limit:=2700
)

This will get you an array of 2700 points.

  1. Now collect these points (to get geometries instead of an array) and use the function convex_hull( ) - all together it looks like this, thus you could skip step 3 and only apply this expression (with Geometry generator or Geometry by expression, see here for details):
convex_hull( 
     collect_geometries (
        overlay_nearest (
            '3000_points',
            $geometry,
            limit:=2700
        )
    )
)

Screenshot: the yellow point (red arrow) is the center point, the red polygon the convex hull covering 2700 of 3000 points: enter image description here

Onion peeling set of minimun-area polygons:

You can now use array_foreach () to create several polygons, containing an ascending percentage of points of the layer. Let's say you want polygons from 1% to 100% in 1% steps. As first argument, create an array like (1,2,3,...., 99, 100) - the easiest way to do so is the function generate_series (1, 100, 1). For each element, you now create a polygon, where the limit argument of overlay_nearest () is the whole number of features divided by 100 and muliplied by the current value of the array - thus for 1: 3000/1001=30 (1%), for 2: 3000/1001=60 (2%) etc.

To convert the resulting array to geometries, again use collect_geometries (). All together, the expression looks like this:

collect_geometries (
    array_foreach (
        generate_series (1, 100, 1),
        convex_hull( 
            collect_geometries (
                overlay_nearest (
                    '3000_points',
                    $geometry,
                    limit:=3000/100*@element
                )
            )
        )
    )
)

Screenshot: a polygon for each integer percentage from 1 to 100: enter image description here

As you can see, the smallest polygon is probably not the smallest possible containing 1% of all points, thus 30 points. There are local clusters where points are very dense, so to include 30 points, there you could create smaller hulls. To realise that, you have again to repeat step 1 from above, but this time with a limit of 30.

So the center point is not the same for each percentage. To illustrate this, I run the whole process for 30 points (1%, red polygon), 100 points (3.33%, blue polygon) and 1000 points (33.33%, yellow polygon) - see the screenshot:

enter image description here

Edit to answer a question from a comment:

Why do I start in the middle? This in fact is not a decision I made, but is inherent to the problem. If you want the convex hull with the smallest possible area, then automatically points near the middle are in a better position. From a corner, to cover 90% of the points, you must reach almost to the other end of the extent of the point layer, thus there is a huge distance. From the corner, you only have neighboring points in one direction. From the middle, you have neighboring points in all direction. Consider drawing circles (buffers) around each point - if you have an even or random distribution of the points, points in the middle have a higher chance to cover more neighboring points with the same radius than points on the corner.

Consider the following example with 50 random points. For each point, I calculated the accumulated distance to the 45 nearest points as explained above. The smallest value I get for the red dot, the highest one for the blue one. Compare their convex hull, covering each time 45 points (90%). Area measurement confirms: the area of the blue polygon is 9% larger than the one of the red one, even though both cover 45 points:

enter image description here

Starting from the corner would only be an option if you have an extreme high density concentration (clustering) of almost all points near a corner with only few points loosely distributed over the rest of the extent. But even in this case, the solution described above will automatically return the "best" center point (of starting point would be a better name for it) from where to start creating minimum convex hull. Also in this case, it center/starting point will be somewhere in the middle of the clustered points, not at the corner - for the same reason explained above.

Accuracy

This method is not 100% accurate. But it is a heuristic approach, a quick and efficient way to get close to the "perfect" solution, getting close to a minimum convex hull with reasonable effort. As @wingnut mentioned in a comment, there probably is no solution to do it deterministically in reasonable time.

To illustrate this with the 50 points from above: For each of the 50 points, I created the convex hull covering the nearest 45 points (90%) with geometry by expression and than calculated the area of the resulting polygons.

The smallest polygon has an area of 3.672, the second smallest 3.679, third smallest 3.679, fourth smallest 3.679, fifth smallest 3.722: this is the red polygon from above. So the solution did not find the smallest, but only the fifth smallest with an area that is 1.3% bigger than the "perfect" solution. The value for the highest area is 4.25 - thus 15.7% larger than the smallest one (the "perfect" solution). So the automatical solution in this case out of 50 possible convex hulls (polygons) found the 5th best solution with a deviation of 1.3% from the best solution.

The blue polygon from above, by the way, is the 11th largest of the convex hulls with an area of 4.060. So it's not the worst, but only the 11th worst (from 50) possibilities, ca. 4.5% smaller than the largest polygon.

5
  • 1
    You are starting in the middle. Why not start in the corner? What is special about the middle.
    – wingnut
    Apr 26 at 8:08
  • Added to my answer a section at the bottom to reply to this question. It's not me to start in the middle, the "ideal" solution, creating a minimum area for the convex hull, tells me to do so. You can test it with own random points and applying the solution above: create convex hulls for the point with the smallest and largest value as calculated in step 1: the one with the smallest value will also rusult in a smaller area in will be located somewhere in the middle. The one with the highest value will almost certainly be located somewhere towards the corner and result in a higher area.
    – Babel
    Apr 26 at 8:52
  • 1
    I like your solution. It is informative. It works particularly well on evenly distributed data.
    – wingnut
    Apr 26 at 9:16
  • Exactly - creating random points results in more or less evenly distributed points. Would be interesting to see the effects for un-evenly distributed data.
    – Babel
    Apr 26 at 9:18
  • 2
    I deal with one particular np-complete problem (an asymmetric TSP), and have had to resort to genetic and reinforcement-learning algorithms to get approximate (but reasonable) solutions in reasonable time. One advantage of a reasonable-looking approximate solution is that it can be considered optimal until and unless someone finds a better one!
    – wingnut
    Apr 26 at 9:22

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.