I've a calculated a triangulation of a point set grouped in certain distance (buffer). For areal classification purposes, I need to calculate:
(a) the points in the unified buffer, which is easy
(b) and need to remove points from the mesh if they are isolated, which is difficult for me with postgis board utilities.
I try to find a function which calculates the edge length for neighbored points (delauney triangulation) and remove points/triangles from the triangulation with large edges.
Details:
Assuming I've a group of points in a plane. The points should be grouped by a set of rules.
1. Collection step Each point with an overlapping buffer of 75m is belonging to a group.
-- 1. Add a geometry buffer for each point, here 75 m
SELECT AddGeometryColumn ('public','stones','b75',32632,'POLYGON',2);
UPDATE stones set b75=ST_SETSRID(ST_BUFFER(geom,75),32632);
-- 2. Add a column rid for a region identifier
ALTER TABLE stones ADD COLUMN rid INT;
-- 3. Create a region table for the overlapping
-- objects of the points
SELECT * INTO stones_bu75 FROM
(SELECT
(A.D).path[1] AS id, -- region identifier
(A.D).geom AS geom -- unified point buffers
FROM
(SELECT st_dump(st_union(bu75)) AS D FROM stones) AS A ) AS S;
-- 4. update the table of points (stones) and mark/collect each
-- stone into the region group
UPDATE stones
SET rid = C.RID FROM -- set the region identifier
( SELECT
S.id AS SID ,
R.id AS RID
FROM
stones AS S, -- get the stone position
stones_bu75 AS R -- get the region
WHERE
-- check if the stone position is in the region
ST_CONTAINS(st_setsrid(B.geom,32632),
st_setsrid(S.geom,32632))
) AS C
WHERE
id = C.SID;
As a result I've stone positions grouped by a region id
and can count them (or sort sparse regions out).
2. Compacting step: In a second step I want to sort out points that are isolated, which means a group of points is valid if there are more than 20 points group within a edge length of < 50m to the next neighbor.
So I need to calculate the triangulation within the groups and filter out the large edges:
-- 1. Triangulation of the points for each group
SELECT AddGeometryColumn ('public','stones_bu75','t75',32632,'TIN',2);
UPDATE stones_bu75 SET t75= T.tri FROM
(SELECT
rid AS tid,
ST_SETSRID(ST_DelaunayTriangles(ST_COLLECT(geom),1, 2),32632) as tri
FROM
(SELECT rid, array_agg(geom) as geom FROM stones GROUP BY rid) AS G
) AS T
WHERE id = tid;
-- 2. Edge length calculation and filtering
geo=# SELECT id, p, ST_length( ST_MakeLine(sp,ep) )
FROM
-- extract the endpoints for every 2-point line segment for each linestring
(SELECT id, p,
ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) as sp,
ST_PointN(geom, generate_series(2, ST_NPoints(geom) )) as ep
FROM
-- extract the individual linestrings
(SELECT id, (G.D).path[1] as p, (ST_Dump(ST_Boundary((G.D).geom))).geom
FROM
(SELECT id, st_dump(t75) AS D from stones_bu75) G
) AS linestrings
) AS segments WHERE ST_length( ST_MakeLine(sp,ep)) < 50;
Question:
Unfortunately, there is no back references to the id of the points spanning the triangulation (only coordinates). Are there possibilities:
(a) to have a logical access (id) to the points and neighbors of each point from the triangulation set,
(b) how do I remove isolated points (points with no connection to the mesh) from a calculated triangulation.
SELECT *, ST_ClusterDBSCAN(geom, 50, 20) OVER() AS clst_id FROM stones;
is what you are after; creating clusters of minimum of 20 points that are no more than 50 meters apart? Note that your points need to be projected into a meter based projection to use 50 (meter) aseps
distance. Each row will get aclst_id
with either a number indicating the same cluster, orNULL
if no cluster was found (isolated)