6

I have some PostGIS SQL that generates uniformly distributed point cluster centroids within polygons (loosely based on this tutorial). The number of point cluster centroids generated is based on the polygon area. The purpose of this is to create uniformly distributed sample points by area within forest stand polygons.

Here is an example of the output:

enter image description here

And here is the SQL (including sample data) to generate this:

-- Make up some data
CREATE TABLE polys(poly_id, geom) AS (
        VALUES  (1, 'POLYGON((1 1, 1 5, 4 5, 4 4, 2 4, 2 2, 4 2, 4 1, 1 1))'::GEOMETRY),
                (2, 'POLYGON((6 6, 6 10, 8 10, 9 7, 8 6, 6 6))'::GEOMETRY)
    );

-- Create point clusters within each polygon
CREATE TABLE pnt_clusters AS
  SELECT  polys.poly_id,
      CASE
          WHEN ST_Area(polys.geom) >9 THEN ST_ClusterKMeans(pts.geom, 8) OVER(PARTITION BY polys.poly_id) 
          
          ELSE ST_ClusterKMeans(pts.geom, 2) OVER(PARTITION BY polys.poly_id) 
         
      END AS cluster_id,
        pts.geom
  FROM    polys,
          LATERAL ST_Dump(ST_GeneratePoints(polys.geom, 1000, 1)) AS pts
;
 
-- Create centroids from point clusters, making sure all centroids are within a polygon
CREATE TABLE centroids AS
  SELECT cluster_id, ST_PointOnSurface(ST_collect(geom)) AS geom
  FROM pnt_clusters
  GROUP BY poly_id, cluster_id;

The "C" shaped polygon on the left contains centroids that are too close to the polygon boundary for my needs. Based on feedback, there seem to be two approaches to handle this situation:

  1. Generate a negative buffer of the polygons and generate the centroids within those negative buffer areas.

  2. Snap the centroids that are within 0.4 units of the original polygon boundary to the closest location on a negative buffer boundary

Option 1 will not work because the density of the point cluster centroids will be too biased toward the center of the polygons. For example, if I generated point cluster centroids within the green negative buffer region, you can imagine the points would be biased too far toward the center of the original forest stand (blue) polygons:

enter image description here

How can I snap the point cluster centroids (red points) that are within 0.4 units of their polygon boundary (blue polygons) to the closest location on the negative buffer polygon boundary (green polygon)?

6
  • 8
    Instead of calculating their locations using the polygon boundaries you could calculate them using the polygon boundary of a buffer inside the original polygons. – PolyGeo Oct 27 '20 at 3:44
  • @PolyGeo I would be concerned that the points may become too concentrated toward the center of the polygon if I used that approach. I was considering using a negative buffer and snapping points close to the border to that negative buffer. I recall hearing a PostGIS pro saying that most operations like this should not be using buffer operations, but rather some sort of distance calculation for the sake of efficiency. – Aaron Oct 27 '20 at 4:04
  • 2
    Go with a negative buffer and produce the sample within it; snapping points only reduces entropy of your samples and increases relative density just as well. If a minimal distance to your borders is required, it is not only unavoidable, but preferable to increase relative density. That, or reduce sample size, e.g. 1000 * ST_Area(<negative_buffer>) / ST_Area(<original_polygon>). – geozelot Oct 27 '20 at 9:30
  • 1
    If I understood the question correctly, you can snap your points on the nearest buffer position with the function ST_ClosestPoint – Atm Nov 10 '20 at 12:07
  • 1
    And in this case it is necessary to think how to make the ST_ClosestPoint function does not create additional points, i.e. 1 centroid from the "centroids" table must be bound by 1 point and not by 2 points to the boundary of the negative buffer... – Cyril Mikhalchenko Nov 10 '20 at 12:57
5
+100

I publish my answer as a workaround to get the expected result for your situation, in the hope that someone will clarify the behavior of ST_ClosestPoint.

So run the script in the form of CTE:

CREATE TABLE polys_pts AS
WITH  
polys(poly_id, geom) AS (VALUES  (1, 'POLYGON((1 1, 1 5, 4 5, 4 4, 2 4, 2 2, 4 2, 4 1, 1 1))'::GEOMETRY),
                (2, 'POLYGON((6 6, 6 10, 8 10, 9 7, 8 6, 6 6))'::GEOMETRY)),
pnt_clusters AS (SELECT  polys.poly_id,
      CASE
          WHEN ST_Area(polys.geom)>9 THEN ST_ClusterKMeans(pts.geom, 8) OVER(PARTITION BY polys.poly_id) 
          ELSE ST_ClusterKMeans(pts.geom, 2) OVER(PARTITION BY polys.poly_id) 
      END AS cluster_id, pts.geom FROM polys,
          LATERAL ST_Dump(ST_GeneratePoints(polys.geom, 1000, 1)) AS pts),
centroids AS (SELECT cluster_id, ST_PointOnSurface(ST_collect(geom)) AS geom FROM pnt_clusters GROUP BY poly_id, cluster_id),
neg_buffer AS (SELECT poly_id, ST_ExteriorRing(ST_Buffer(geom, -0.4, 'endcap=flat join=round')) geom FROM polys GROUP BY poly_id, polys.geom),
pnt_clusters_new AS (SELECT DISTINCT ST_ClosestPoint(a.geom, b.geom) AS geom FROM neg_buffer a, centroids b),
node_pts AS (SELECT ST_StartPoint(geom) geom FROM neg_buffer),
snap_pts AS (SELECT b.cluster_id, a.geom FROM pnt_clusters_new a JOIN centroids b ON ST_DWithin(a.geom, b.geom, 0.4))
SELECT  a.cluster_id, (a.geom) geom FROM snap_pts a WHERE NOT EXISTS (SELECT 1 FROM node_pts b WHERE ST_Intersects(a.geom, b.geom))

And check the result.


I hope I have understood your qualifying condition correctly and as a result run the specified script:

CREATE TABLE polys_pts AS
WITH  
polys(poly_id, geom) AS (VALUES  (1, 'POLYGON((1 1, 1 5, 4 5, 4 4, 2 4, 2 2, 4 2, 4 1, 1 1))'::GEOMETRY),
                (2, 'POLYGON((6 6, 6 10, 8 10, 9 7, 8 6, 6 6))'::GEOMETRY)),
pnt_clusters AS (SELECT  polys.poly_id,
      CASE
          WHEN ST_Area(polys.geom)>9 THEN ST_ClusterKMeans(pts.geom, 8) OVER(PARTITION BY polys.poly_id) 
          ELSE ST_ClusterKMeans(pts.geom, 2) OVER(PARTITION BY polys.poly_id) 
      END AS cluster_id, pts.geom FROM polys,
          LATERAL ST_Dump(ST_GeneratePoints(polys.geom, 1000, 1)) AS pts),
centroids AS (SELECT cluster_id, ST_PointOnSurface(ST_collect(geom)) AS geom FROM pnt_clusters GROUP BY poly_id, cluster_id),
neg_buffer AS (SELECT poly_id, (ST_Buffer(geom, -0.4, 'endcap=flat join=round')) geom FROM polys GROUP BY poly_id, polys.geom),
neg_buffer_pts_out AS (SELECT a.cluster_id, (a.geom) geom FROM centroids a WHERE NOT EXISTS (SELECT 1 FROM neg_buffer b WHERE ST_Intersects(a.geom, b.geom))),
neg_buffer_pts_in AS (SELECT a.cluster_id, (a.geom) geom FROM centroids a WHERE EXISTS (SELECT 1 FROM neg_buffer b WHERE ST_Intersects(a.geom, b.geom))),
snap_pts_clusters_in AS (SELECT DISTINCT ST_ClosestPoint(ST_ExteriorRing(a.geom), b.geom) AS geom FROM neg_buffer a, neg_buffer_pts_in b),
node_pts AS (SELECT ST_StartPoint(ST_ExteriorRing(geom)) geom FROM neg_buffer),
snap_pts AS (SELECT b.cluster_id, a.geom FROM snap_pts_clusters_in a JOIN centroids b ON ST_DWithin(a.geom, b.geom, 0.4))
SELECT  a.cluster_id, (a.geom) geom FROM snap_pts a WHERE NOT EXISTS (SELECT 1 FROM node_pts b WHERE ST_Intersects(a.geom, b.geom))
UNION SELECT c.cluster_id, (c.geom) geom FROM neg_buffer_pts_out c ORDER BY cluster_id;

Check the result and refine its behavior if necessary...

As you can see from the script, I added a part of the code that separated the points by the hit and miss condition (buffer zone), and the logic remained the same.

However, as a result, my new question appeared about the strange behavior of the ST_ClosestPoint() function.

2
  • This is a nice approach! However, I do notice that it snapped all of the point cluster centroids to the closest negative buffer location, rather than only the points that are within 0.4 units of the original polygon boundary (i.e. the two points in the "C" shaped polygon). Any thoughts on how to handle that? – Aaron Nov 11 '20 at 2:36
  • 1
    Looking at this and similar geoinstruments, I exclaim: "What interesting results can a collective creative work lead to"! Such behavior is the future of geoinformation development (despite the programming language)...This situation is associated with a long distance team race, the first ran "geozelot", then he passed the baton "Aaron", then I got the baton and ran to the finish line, a funny story has happened...although it is necessary to look and wait, maybe someone is already running next to me, and the finish is still far away :-)...a joke. – Cyril Mikhalchenko Nov 11 '20 at 12:30
4

There will be bias, no matter what you do.

ST_PointOnSurface operates with a deterministic algorithm based on the areal shape; in combination with your cluster based centroids, this has the effect that, when the area of the ST_Area(ST_Buffer(geom, -0.4)) is less than 50% of the ST_Area(geom), you may get POS outside of that threshold. However, since the general extent and shape of the areas that make up each cluster stays similar despite scaling the initial polygons, the general distribution of centroid would also stay similar, even if you scale input area for ST_GeneratePoints. For centroids that would fall outside of the actual area (like with the C shape), the POS would be placed on the boundary if the boundary is more than half way away from the centroid.

In other words, you should get equally uniformly distributed centroids, always within the threshold distance, when creating random points in a negatively buffered polygon only if the buffered area would be less than half the size of the actual area, and else use the actual area:

SELECT  poly_id,
        cluster_id,
        ST_PointOnSurface(ST_Collect(geom)) AS geom
FROM    (
SELECT  polys.poly_id,
        CASE WHEN ST_Area(polys.geom) > 9
            THEN ST_ClusterKMeans(pts.geom, 8) OVER(PARTITION BY polys.poly_id) 
            ELSE ST_ClusterKMeans(pts.geom, 2) OVER(PARTITION BY polys.poly_id) 
        END AS cluster_id,
        pts.geom
FROM    polys,
        LATERAL ST_Dump(ST_GeneratePoints(CASE ST_Area(ST_Buffer(polys.geom, -0.4)) / ST_Area(polys.geom) > 0.5 WHEN TRUE THEN geom ELSE ST_Buffer(polys.geom, -0.4) END, 1000, 1)) AS pts
) q
GROUP BY
        1, 2
;

enter image description here

[ lighter gray: initial shapes; darker gray: negative buffer (-0.4); red: no buffer; yellow: ratio dependent buffer; orange: points fall on the same spot]

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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