11

I've to create dissolved buffers from multi-point input features. In the example below, the input table contains 4 features. Feature #2 consists of two point geometries. After creating a buffer, I get 4 polygon geometries:

enter image description here

Is there a way to group the result? The buffers of the points #1 and #2 are dissolved and should be a single multi-polygon feature (a).

What I've done so far:

-- collect all buffers to a single multi-polygon feature
-- dissolve overlapping polygon geometries
CREATE TABLE public.pg_multibuffer AS SELECT
    row_number() over() AS gid,
    sub_qry.*
FROM (SELECT
    ST_Union(ST_Buffer(geom, 1000, 8))::geometry(MultiPolygon, /*SRID*/) AS geom
FROM
public.multipoints)
AS sub_qry;

EDIT:

-- create sample geometries

CREATE TABLE public.multipoints (
gid serial NOT NULL,
geom geometry(MultiPoint, 31256),
CONSTRAINT multipoints_pkey PRIMARY KEY (gid)
);

CREATE INDEX sidx_multipoints_geom
ON public.multipoints
USING gist
(geom);

INSERT INTO public.multipoints (gid, geom) VALUES
(1, ST_SetSRID(ST_GeomFromText('MultiPoint(12370 361685)'), 31256)),
(2, ST_SetSRID(ST_GeomFromText('MultiPoint(13520 360880, 19325 364350)'), 31256)),
(3, ST_SetSRID(ST_GeomFromText('MultiPoint(11785 367775)'), 31256)),
(4, ST_SetSRID(ST_GeomFromText('MultiPoint(19525 356305)'), 31256));
5
  • You use subqueries too much. This elimininates your ability to GROUP BY on the attribute on which you want to cluster.
    – Vince
    Mar 26, 2017 at 12:27
  • So you need to do a spatial union and then also a union based on feature number, which is why you expect 3 multipolygons from the diagram above. I suspect this will require a two-step process, but just wanted to be clear on the question, before offering an answer. Mar 26, 2017 at 20:58
  • Yes, I want to union the buffer polygons and collect the result based on the number of the input features. Mar 26, 2017 at 21:51
  • Any update on this? I would like to know if this works for you, as as far as I can see, I have answered the question. Apr 19, 2017 at 16:31
  • Sorry for the late reply, I haven't been online for a couple of days. Apr 21, 2017 at 21:25

3 Answers 3

8

Starting with some random points, in an attept to imitate those in the OP's image, where the first two spatially intersect, then the 2nd and 3rd have the same attribute id (2), with a couple of other points that neither spatially intersect nor have the same attribute, the following query produces 3 clusters:

WITH 
  temp (id, geom) AS 
     (VALUES (1, ST_Buffer(ST_Makepoint(0, 0), 2)),
        (2, ST_Buffer(ST_MakePoint(-0.7,0.5), 2)),
        (2, ST_Buffer(ST_MakePoint(10, 10), 2)), 
        (3, ST_Buffer(ST_MakePoint(-2, 12), 2)), 
        (4, ST_Buffer(ST_MakePoint(5, -6), 2))),
 unions(geoms) AS 
      (SELECT ST_Union(geom) FROM temp GROUP BY id),
 clusters(geoms) AS 
      (SELECT ST_CollectionExtract(unnest(ST_ClusterIntersecting(geoms)), 3) 
         FROM unions),
 multis(id, geoms) AS 
      (SELECT row_number() over() as id, geoms FROM clusters)
 SELECT ST_UNION(d.geom) FROM 
      (SELECT id, (ST_DUMP(geoms)).geom FROM multis) d GROUP BY id;

There are several steps here:

  1. use ST_Union, grouping by id, to first group by attribute
  2. use ST_ClusterIntersecting to combine those from same group that intersect spatially
  3. add an id to each of the clusters (table multis) -- trying to do this directly in the ClusterIntersecting leads to all geometries getting an id of 1
  4. Union the dumped geometries from step 2, grouping by the id from step 3 -- this is the dissolve part. This causes the two overlapping polygons in your cluster A, to be joined together, rather than being overlapping, as they are at the end of step 2.

Rather long, but it works (and, I am sure there is a shorter way).

Using the WKT tool in QGIS, (and discovering how awful I am with the editing tools) produces clusters like the following, where you can see the cluster your- labelled as a, is all together -- ie, one colour.

enter image description here

If you put an ST_AsText round the final, ST_UNION(d.geom), then you can see the results directly.

EDIT following more information in the comments: As you are starting with points you will need to incorporate the buffer into my original solution -- which I put in the temp CTE at the start to mimic your diagram. It would be easier to add the buffer in the unions CTE, so you can do all the geometries at once. So, using a buffer distance of 1000, as an example, the following now returns 3 clusters, as expected.

WITH temp(id, geom) AS 
  (VALUES 
      (1, ST_SetSRID(ST_GeomFromText('MultiPoint(12370 361685)'), 31256)),   
      (2, ST_SetSRID(ST_GeomFromText('MultiPoint(13520 360880, 19325 364350)'), 31256)),                                                
      (3, ST_SetSRID(ST_GeomFromText('MultiPoint(11785 367775)'), 31256)),
      (4, ST_SetSRID(ST_GeomFromText('MultiPoint(19525 356305)'), 31256))
),                                              
unions(geoms) AS 
  (SELECT st_buffer(ST_Union(geom), 1000) FROM temp GROUP BY id),
clusters(geoms) AS 
  (SELECT ST_CollectionExtract(unnest(ST_ClusterIntersecting(geoms)), 3) 
     FROM unions),
multis(id, geoms) AS 
  (SELECT row_number() over() as id, geoms FROM clusters)
SELECT id, ST_UNION(d.geom) FROM 
  (SELECT id, (ST_DUMP(geoms)).geom FROM multis) d GROUP BY id;
4
  • Sorry it took so long to get back to you. I have some trouble to visualize the buffer geometries in QGIS. I have tried to modify your query using ST_SetSRID, ST_Multi and ::geometry(Multipolygon, /*SRID*/), but at the moment it‘s not working. Apr 3, 2017 at 12:45
  • OK, if you can post your code, and even better some data, I might be able to help. Apr 3, 2017 at 12:52
  • I've added some SQL to create sample points. Apr 3, 2017 at 14:05
  • Bit tied up today, I'll revert as soon as I can. Will have to work the multipoint into the query too. Apr 3, 2017 at 14:59
3

One way to do this is to ST_Union all of the buffers together, ST_Dump the result to get the components of the resulting polygon, and join with ST_Intersects back to the input points to found out how many/which points made up each cluster.

This can be done without requiring a join by grouping the points together before calling ST_Buffer. For two points to be located within the same dissolved buffer, they must be reachable by hops between points of a distance less than eps. This is just a minimum-linkage clustering problem, which can be solved using ST_ClusterDBSCAN:

SELECT
  cluster_id,
  ST_Union(ST_Buffer(geom, 1000)) AS geom,
  count(*)                        AS num_points,
  array_agg(point_id)             AS point_ids
FROM (
  SELECT
    point_id,
    ST_ClusterDBSCAN(geom, eps := 2000, minpoints := 1) OVER() AS cluster_id ,
    geom
  FROM points) sq
 GROUP BY cluster_id;

Note that this won't produce exactly the same result as the buffer-first method, because PostGIS buffers aren't perfect circles and two points 1000m apart may not be connected by two 500m buffers.

2
  • Seems we had a similar idea. I haven't tested yours, but I'm sure it works, and more cleanly than mine. Mar 27, 2017 at 14:00
  • It seems that PostGIS 2.2.1 doesn't support ST_ClusterDBSCAN. I‘ve installed PostGIS 2.3.2, but new postgis extensions in pgAdmin are still version 2.2.1. Apr 2, 2017 at 20:04
1

According to this answer you want to do ST_DUMP within your subquery.

Something like this:

-- collect all buffers to a single multi-polygon feature
-- dissolve overlapping polygon geometries
CREATE TABLE public.pg_multibuffer AS SELECT
    row_number() over() AS gid,
    sub_qry.*
FROM (SELECT
    ST_Dump(ST_Union(ST_Buffer(geom, 1000, 8))::geometry(MultiPolygon, /*SRID*/)) AS geom
FROM
public.multipoints)
AS sub_qry;

The reason why is that ST_UNION returns a dissolved multipolygon of all the features, and ST_DUMP breaks this out into the individual polygon features (which were dissolved).

3
  • 1
    This won't actually work, because any attributes which would be needed to cluster the desired multipart polygon will have been lost.
    – Vince
    Mar 26, 2017 at 14:04
  • I've tried ST_Multi((ST_Dump(ST_Union(ST_Buffer(geom, 1000, 8)))).geom)::geometry(MultiPolygon, /*SRID*/) AS geom, but this creates 4 features instead of 3. Mar 26, 2017 at 14:43
  • Oh, right, you want to group by the number? You'll need to GROUP_BY before you ST_UNION.
    – Alex Leith
    Mar 27, 2017 at 1:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.