# How to group polygons based on spatial proximity and field value?

I have a shapefile consisting of the day of the year on which the temperature in a lake exceeds 15 degrees Celsius. I can see clear patterns in this (here seen as a raster):

I want to group the pixels adjacent to each other if their difference is no larger than 5 days. This means, I both have to select a spatial join but also group via a field value - in this case day of the year `doy`. Eventually, it should look a bit like this (conceptual) with the colors indicating different groups, and the numbers representing day of the year:

I have tried to use the `aggregate` tool in QGIS with the expression `floor((doy - 66)/5)`, but here I create 5 groups. This means I will have pixel values that are in close proximity - spatially and in value - but be split into two different groups as it is based on 5 static groups.

I would like to add the pixels together in a group by spatial proximity and a "running" statement being something like `if doy of pixel A < 5 doy pixel B then group`.

Is it possible in QGIS?

• So one group could be 100-104-108-112-116-120-124? Or should the new pixel be compared to the average of the existing pixels in that group?
– Bera
Commented Jun 29, 2021 at 12:21
• @BERA - that was my original idea here: gis.stackexchange.com/a/402492/88814 - but it seems not to be what OP is looking for, see comments to my solution. Commented Jun 29, 2021 at 12:39

The solution: principles

Automatically draw a connecting each cell to each neighboring cell with a `doy`-value in the range of +/-5 using QGIS expressions with Geometry by expression. Buffer and dissolve the line to get a small, line-like polygons that connect all the cells that should be grouped together. Now assign to each grid-cell the id of the feature of the line-polygon that covers it. Then group the cells together using `Aggregate` by this id.

The workflow

You can use QGIS expression using `overlay_nearest()` (available since QGIS 3.16) to get for each cell the id of all neighboring cells that have a `doy` value that is in the range of +/-5 to the `doy`-value of the current cell. Use these values to create a line connecting the centroid of each cell to each neighboring cell's centroid that fulfills that condition. See below for the expression that creates these lines.

When you have these lines (see screenshot 1) use `Select by expression` with this expression `\$length=0` to select all lines with a length of 0 and delete them.

Then buffer the remaining lines with a very small distance (clearly smaller than the cell size) and dissolve the buffer. Then convert multi- to singleparts. Now you have a small, line-like polygon connecting all the cells you want to group together. Creat a new field named `agg_id` with unique values using field calculator with the expression `\$id`.

Now `Join attributes by location` to the grid to assign the `agg_id` from the buffer-polygon to the grid. Then use `Aggregate` to get the result you're looking for (screenshot 2). You're done.

Screenshot 1, exemplifying the connection of neighboring cells if the difference of their `doy`-value is smaller than 5, using a grid with random values:

Screenshot 2: the result of this solution, based on the same data, showing the grouped together (aggregated) cells in a red frame:

Details: the expression

1. From the grid, extract centroids.

2. On this centroid layer, use the expression from below (after the screenshot) with `Geometry by expression` to create the lines that connect the cells you want to group together. Replace `centroid` with the name of your centroid-layer.

On this screenshot, you see the expression in action (here used with geometry generator). Bold label in black: `id` and `doy` value of each cell, in smaller red script the id of those neighboring cells that should be grouped to the current cell:

``````collect_geometries (
array_remove_all(
with_variable (
'line',
'make_line (\$geometry,geometry (get_feature_by_id (''centroid'',@element)))',
array_foreach (
array_remove_all(
array_foreach (
overlay_nearest(
'centroid',
id,
limit:=5
),
if (
array_contains(
array_remove_all(
array_foreach (
overlay_nearest(
'centroid',
doy,
limit:=5
),
if (
abs(doy-@element)<=5,
@element,
'del'
)
),
'del'
),
attribute(
get_feature_by_id (
'centroid',
@element
),
'doy'
)
)=true,
@element,
'del'
)
),
'del'
),
if (
degrees (
azimuth (
start_point (eval(@line)), end_point (eval(@line))
)
) % 90 > 5,
'del',
eval(@line)
)
)
),
'del'
)
)

``````
• It is a classic task for networkx module, which has connected components function. 5 lines of code, when links between nodes established. Commented Jun 29, 2021 at 21:31
• OK, would be nice to see it added as another solution. Commented Jun 29, 2021 at 21:32
• I can post arcmap solution, I don't think it's OK for qgis tag Commented Jun 29, 2021 at 21:45
• Maybe you can explain the conceptual principles - this should be possible to adapt in QGIS? Commented Jun 29, 2021 at 21:54
• I like this solution - it works with pixels that are aligned and the same size. However, with dissolved pixels, like mine, some centroids are not fully aligned and, therefore, the line drawn cannot connect to those. So it did not fully resolve my problem, but may solve other people's problems. Commented Jul 1, 2021 at 11:21

Solution for ArcMap. Picture shows original polygons, called NODES, and values they store:

Table shows NODES table, highlighted column stores their names:

Table, called EDGES shows numerically connections between polygons (I used Polygon Neighbors tool to compute it):

Highlighted columns store NODES' names. At this stage we can simply remove records where

``````abs(src_VALUE-nbrVALUE) > 5
``````

Script below assumes we've done that.

``````import arcpy
import networkx as nx

G=nx.Graph() #  create undirected graph
EDGES = arcpy.da.TableToNumPyArray("EDGES", ("src_FID","nbr_FID")) #  read names of nodes connected to each other

dictFeatures = {} # dictionary to populate by group number
for m,group in enumerate(nx.connected_components(G)):
for node in group:dictFeatures[node]=m

with arcpy.da.UpdateCursor("NODES",("FID","PART")) as cursor: # Transfering results to NODES
for fid,prt in cursor:
prt=dictFeatures.get(fid,-1)
cursor.updateRow((fid,prt))
``````

NODES colored by their group number:

Note: script populates result field (called PART in our example) by -1 for 'island' nodes.