Spatial Clustering - associate a cluster attribute (id) with geometry that is part of the cluster
I am having some problems to relate the clustered set of geometries to their own propriety.
Data
I have a table with a set of geometries,
buildings {
gid integer,
geom geometry(Multipoligon,4326)
}
And I ran the ST_ClusterWithin function with a specific threshold over the "buildings" table. From this cluster analysis I got a table called "clusters"
clusters {
cid Integer,
geom geometry(GeometryCollection,4326)
}
Question
I would like to extract into a table all the geometry with its own corresponding cluster information.
clustered_building {
gid Integer
cid Integer
geom geometry(Multipoligon,4326)
}
gid | cid | geom |
-----+------------+-----------------------+
1 | 1 | multypoligon(...) |
2 | 1 | multypoligon(...) |
3 | 1 | multypoligon(...) |
4 | 2 | multypoligon(...) |
5 | 3 | multypoligon(...) |
6 | 3 | multypoligon(...) |
What I did (but doesn't work)
I tried to use two ST_GeometryN / ST_NumGeometries functions analyze each MultyGeometry and extract the cluster information with this query, which is derived from one of the standard ST_Geometry man page examples.
INSERT INTO clustered_building (cid, c_item , geom)
SELECT sel.cid, n, ST_GeometryN(sel.geom, n) as singlegeom
FROM ( SELECT cid, geom, ST_NumGeometries(geom) as num
FROM clusters") AS sel
CROSS JOIN generate_series(1,sel.num) n
WHERE n <= ST_NumGeometries(sel.geom);
It takes a few seconds in the request if I am forced to use a series of 10.
CROSS JOIN generate_series(1,10)
But it got stuck when I ask to generate a series according to the number of elements in each GeometryCollection. Also, this query does not allow linking a single geometry with its own functions in the build table because I am losing the "gid"
Can someone please help me, thanks
Stefano
source to share
I have no data, but using some dummy values where ids 1, 2 and 3 overlap and 4 and 5, you can do something like the following:
WITH
temp (id, geom) AS
(VALUES (1, ST_Buffer(ST_Makepoint(0, 0), 2)),
(2, ST_Buffer(ST_MakePoint(1, 1), 2)),
(3, ST_Buffer(ST_MakePoint(2, 2), 2)),
(4, ST_Buffer(ST_MakePoint(9, 9), 2)),
(5, ST_Buffer(ST_MakePoint(10, 10), 2))),
clusters(geom) as
(SELECT
ST_Makevalid(
ST_CollectionExtract(
unnest(ST_ClusterIntersecting(geom)), 3))
FROM temp
)
SELECT array_agg(temp.id), cl.geom
FROM clusters cl, temp
WHERE ST_Intersects(cl.geom, temp.geom)
GROUP BY cl.geom;
If you end up with the final cl.geom, which is ST_AsText, you will see something like:
{1,2,3} | MultiPolygon (((2,81905966523328 0.180940334766718,2.66293922460509 -0.111140466039203,2.4142135623731 -0.414213562373094,2.11114046603921 -0.662939224605089,1.81905966523328 -0.819059665233282,1.84775906502257 -0.765366864730179,1.96157056080646 -0.390180644032256,2 0,2 3.08780778723872e-16,2 0,2.39018064403226 0.0384294391935396, 2.76536686473018 0.152240934977427.2.81905966523328 0.180940334766718)) ......
{4,5} | MultiPolygon (((10,8190596652333 8.18094033476672,10.6629392246051 7.8888595339608,10.4142135623731 7.58578643762691,10.1111404660392 7.33706077539491,9.76536686473018 7.15224093497743,9.39018064403226 7.03842943919354,9 7,8.60981935596775 7.03842943919354,8.23463313526982 7.15224093497743,7.8888595339608 7.33706077539491,7.58578643762691 7.5857864376269,7.33706077539491 7.88885953396079,7.15224093497743 8,23463313526982
where you can see IDs 1,2,3, below the first multipolygon and 4.5 others.
The general idea is that you are a cluster of data, and then you traverse the returned clusters with the original data, using array_agg to group the ids together so that the returned multipolygons now contain the original ids. Using ST_CollectionExtract with 3 as the second parameter, combined with disest , which splits the geometry collection returned by ST_ClusterIntersecting back into strings, returns each contiguous cluster as a (polygon) polygon. ST_MakeValid is that sometimes when you intersect geometry with other related geometries, such as the original polygons with your clustered polygons, you get strange rounding effects and GEOS error about unsettled intersections, etc.
I recently answered a similar question on gis.stackexchange that you might find useful.
source to share