8

I have a polygon layer. I would like to select the minimum number of neighbouring polygons to give a total area greater than a fixed value. In reality the polygons are land parcels and I would like to know which parcels I need purchase in addition to create a continuous area of, for example 100 hectares.

Is this possible in QGIS using virtual-layers or by querying the spatialite database in the DB Manager? I need to repeat this process for potentially hundreds of sites, hence the need for an automated process. Is there some other way that could be used to achieve this within QGIS such as the processing modeler?

I had thought of creating a buffer of the appropriate area and then seeing which polygons intersect. However, the resulting area could be too small (not acceptable. Which would mean lots of trial and error for each site until I got it right. An area larger than my target value would be acceptable, however, the excess should minimised.

5
  • The answer is Recursion.. if your original polygon + a isn't big enough use o + a to select b, c, d and then test o + a + b, o + a + c and o + a + d recursively to find the minimum number of polygons that meet the size. Be sure to stop your recursion when size is met for aggregate, no further polygons are selected by the aggregate or another satisfactory conclusion is reached. It might be possible in SpatialLite but I have no experience of that; it could definitely be done in pyqgis but I'm a bit too rusty with that to venture an answer. Feb 17, 2017 at 0:59
  • I figured that would be the way forward, if it was javascript no problem. I have no experience of writing python, maybe it's time get my hands dirty. Feb 17, 2017 at 6:43
  • You could create a custom application with very little code in OpenLayers /Leaflet/Mapbox where you click on a parcel and then, using an intersection function, highlight surrounding parcels until the sum of areas > 100 HA. In some cases you will probably even have different possible combinations of surrounding parcels, because parcels tend to me of similar size. In that case, you could visualize the different combinations of parcels that add up to more than 100 HA. Otherwise, you could experiment in PyQGIS Feb 17, 2017 at 6:52
  • 2
    In graph theory terms, you want to construct the adjacency graph for the polygons with area as vertex weight, and then find the subgraph with the fewest vertexes that has the smallest total vertex weight greater than your threshold. You should probably look at python graph theory packages. You might be able to construct the adjacency graph from existing plugin code - I definitely wrote something to do it for my plugin that coloured polygons such that no two adjacent polygons had the same colour, and I think someone has ported that to the latest QGIS...
    – Spacedman
    Feb 17, 2017 at 8:07
  • look for papers on zone design or redistricting to find good algorithms - usually you'd use simulated annealing or tabu search to find good enough local optima for the problem. google.co.uk/search?q=openshaw+zdes
    – Ian Turton
    Feb 17, 2017 at 9:27

1 Answer 1

2

You can do that using QGIS expressions. The idea behind it is to to use overlay_nearest and to use it to create an array of the area of all neighboring polygons. Sort this array largest to smallest, then check how many elements of the array you need to sum up to reach the threshold - this is the total area you want, minus the area you already have (the initial polygon). Render the geometries of these elements to get a selection.

To implement this, you can use this expression. You can adapt it as follows:

  • line 5 (here: 6000): total area you want to get (adding initial + neighboring polygons)
  • line 8 (here:20): the max. number of eighboring polygons to consider. Reduce the number for faster performance, increase it if you need more neighbors to get the area you want.
  • sixth last line (here: 'area'): name of the attribute that contains the value of the area for each polygon. Here, the area is rounded to 1 digit. If you have another number of digits, adapt this number in the next line: round (@element,1).

Here is the expression:

 if (
     is_selected(),
     collect_geometries (
         array_foreach (
             with_variable(
                 'totsize',
                 6000,
                 with_variable(
                     'limit',
                     20,
                     array_slice (
                         array_sort(
                             overlay_nearest (
                                 @layer, 
                                 $area, 
                                 limit:=@limit
                             ),
                             0
                         ),
                         0,
                         @limit+1-array_sum(
                             array_foreach (
                                 generate_series (0,@limit),
                                 array_sum(
                                     array_slice (
                                         array_sort(
                                             overlay_nearest (
                                                 @layer,
                                                 $area, 
                                                 limit:=@limit
                                             ),
                                             0
                                         ),
                                         0,
                                         @element
                                     )
                                 )>@totsize-$area
                             )
                         )
                     )
                 )
             ),
             geometry (
                 get_feature (
                     @layer,
                     'area',
                     round (@element,1)
                 )
             )
         )
     ),
     ''
 )   

Screenshot: the polygon in the middle (area: 1131.5) is selected and the expression is set to a total area of 600. The 8 largest neighboring polygons are added to reach a total area of 6152.3: 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.