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
Jun 29 '21 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. Jun 29 '21 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. Jun 29 '21 at 21:31
• OK, would be nice to see it added as another solution. Jun 29 '21 at 21:32
• I can post arcmap solution, I don't think it's OK for qgis tag Jun 29 '21 at 21:45
• Maybe you can explain the conceptual principles - this should be possible to adapt in QGIS? Jun 29 '21 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. Jul 1 '21 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.