0

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.

point distribution of stones and boulders

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

enter image description here

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.

2
  • I wonder if 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) as eps distance. Each row will get a clst_id with either a number indicating the same cluster, or NULL if no cluster was found (isolated)
    – geozelot
    Mar 25, 2020 at 14:41
  • I've to check this with a new postgis version, because I'm still running my Debian 9 with Postgis 2.1.
    – huckfinn
    Mar 27, 2020 at 13:35

0

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.

Browse other questions tagged or ask your own question.