3

I am using st_concavehull to create boundaries/shapes. This works, however, the shape is not precisely tracing the coordinates. It is more puffy, as if trying to even out or "clean up" the actual shape.

How can I create a more cleaner/neater boundary? I've tried numbers from -1 to 1 and nothing seems to produce the actual shape I'm expecting; they all have curved edges and "cover" the coordinates instead of tracing them.

Example/My Code:

st_concavehull(
                    st_collect(array[
                        st_transform(st_setsrid(st_makepoint(-97.5660461, 30.4894905), 4326),4269),
st_transform(st_setsrid(st_makepoint(-97.5657216, 30.4902173), 4326),4269),
st_transform(st_setsrid(st_makepoint(-97.5608779, 30.4896142), 4326),4269),
st_transform(st_setsrid(st_makepoint(-97.5605001, 30.491422), 4326),4269),
st_transform(st_setsrid(st_makepoint(-97.5588115, 30.4911697), 4326),4269),
st_transform(st_setsrid(st_makepoint(-97.5588262, 30.4910204), 4326),4269),
st_transform(st_setsrid(st_makepoint(-97.5585742, 30.4909966), 4326),4269),
st_transform(st_setsrid(st_makepoint(-97.5578045, 30.4909263), 4326),4269),
st_transform(st_setsrid(st_makepoint(-97.5574653, 30.4908877), 4326),4269),
st_transform(st_setsrid(st_makepoint(-97.5571534, 30.4908375), 4326),4269),
st_transform(st_setsrid(st_makepoint(-97.5560964, 30.4907427), 4326),4269),
...

                    ]), -0.15)  as polygon
            ) concavehull
        ) );

  • A precise shape around points is bound to (human) interpretation; a set of 4 points in a square plus one in the middle can be precisely traced in at least 4 ways connecting all, and even more when you leave points out. I haven't checked your point set, but unless it's a rare edge case, polygonizing points is bound to ambiguity without a predefined 'order'. – geozelot Feb 7 at 20:53
  • @ThingumaBob is there a way to order my coordinates before hand to produce a more accurate shape? – Megan Feb 7 at 23:17
  • ST_ConcaveHull is generally expected to produce a polygon which does not have all the points in its boundary. But in this case the inability to find a tolerance which produces a finer result is probably a limitation of the implementation. A different Concave Hull algorithm might well allow including all the points presented. – dr_jts Feb 8 at 1:10
  • @Megan Generally, a Polygon is defined by at least one closed LineString, so a 'trace' could be defined by point indices (orientation, i.e. clockwise or not, is important for different formats); then, a (closed) LineString could be generated from points in that order, and get polygonized. – geozelot Feb 8 at 11:46
5

One way to produce a polygon which uses all points in a set is to create a star polygon around the centroid of the points. Whether this produces an acceptable result depends on the positions of the points. But in this case it seems like this will do exactly what is required.

The algorithm is simple:

  • Compute the centroid of the set of points
  • For each point in the set, compute the azimuth (angle) from the centroid to the point
  • Join up the points in order of the azimuth and create a polygon from them

Here's sample code to do this for the data given:

WITH pts(pt) AS (VALUES
(st_transform(st_setsrid(st_makepoint(-97.5660461, 30.4894905), 4326),4269) ),
(st_transform(st_setsrid(st_makepoint(-97.5657216, 30.4902173), 4326),4269) ),
(st_transform(st_setsrid(st_makepoint(-97.5608779, 30.4896142), 4326),4269) ),
(st_transform(st_setsrid(st_makepoint(-97.5605001, 30.491422), 4326),4269) ),
(st_transform(st_setsrid(st_makepoint(-97.5588115, 30.4911697), 4326),4269) ),
(st_transform(st_setsrid(st_makepoint(-97.5588262, 30.4910204), 4326),4269) ),
(st_transform(st_setsrid(st_makepoint(-97.5588262, 30.4910204), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5585742, 30.4909966), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5578045, 30.4909263), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5574653, 30.4908877), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5571534, 30.4908375), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5560964, 30.4907427), 4326),4269))
),
centroid AS (SELECT ST_Centroid( ST_Collect(pt) ) AS centroid FROM pts),
line AS (SELECT ST_MakeLine( pt ORDER BY ST_Azimuth( centroid, pt ) ) AS geom
    FROM pts CROSS JOIN centroid),
poly AS (SELECT ST_MakePolygon( ST_AddPoint( geom, ST_StartPoint( geom ))) AS geom
    FROM line)
SELECT geom FROM poly;

and here is the result:

enter image description here

(This would make a nice function for PostGIS, either built-in or as a community contribution in PL/pgSQL)

CREATE OR REPLACE FUNCTION ST_PointsInStarPolygon(
geom GEOMETRY
)
RETURNS GEOMETRY AS
$BODY$
WITH
pts(geom) AS (SELECT (ST_DumpPoints(geom)).geom geom),
centroid AS (SELECT ST_Centroid(ST_Collect(geom)) AS centroid FROM pts),
line AS (SELECT ST_MakeLine(geom ORDER BY ST_Azimuth(centroid, geom)) AS geom
    FROM pts CROSS JOIN centroid),
poly AS (SELECT ST_MakePolygon(ST_AddPoint(geom, ST_StartPoint(geom))) AS geom
    FROM line)
SELECT geom FROM poly;
$BODY$
LANGUAGE SQL

Run:

SELECT ST_PointsInStarPolygon(geom) geom FROM (SELECT ST_Collect(geom) AS geom FROM pts) a;
| improve this answer | |
  • +1. Then what's the problem? If you don't mind wrapping it in a function, give it a name like ST_PointsInStarPolygon or whatever you like and publish it, it deserves attention. – Cyril Mikhalchenko Feb 8 at 14:42
  • 1
    No problem, just letting someone else have some fun too. If not I'll get there eventually. – dr_jts Feb 8 at 16:32
  • @dr_jts Ah, thank you I am starting to see the shapes better! It seems though that there may be a few straggling polygons who have star shapes sticking out. Is there a way to clean these up. Here is a picture from my results – Megan Feb 11 at 21:08
  • 1
    It would be easier to suggest options if you post sample data that is exhibiting the issues you're seeing. But it looks like there are points quite far inside the polygons, which will of course create "gores". Cleaning this up is going to be art not science. It might be possible to remove narrow reflex vertices. – dr_jts Feb 12 at 22:32
1

Here's another simple, customizable geoinstrument:

Run Postgre/PostGIS SQL in PGAdmin:

WITH pts(geom) AS (VALUES
(st_transform(st_setsrid(st_makepoint(-97.5660461, 30.4894905), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5657216, 30.4902173), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5608779, 30.4896142), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5605001, 30.491422), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5588115, 30.4911697), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5588262, 30.4910204), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5588262, 30.4910204), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5585742, 30.4909966), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5578045, 30.4909263), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5574653, 30.4908877), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5571534, 30.4908375), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5560964, 30.4907427), 4326),4269))), 
tbla AS (SELECT ST_Centroid(ST_Collect(geom)) geom FROM pts),
tblb AS (SELECT geom FROM pts UNION SELECT geom FROM tbla),
tblc AS (SELECT ST_MakeLine(a.geom, b.geom) geom FROM tbla a CROSS JOIN pts b),
tbld AS (SELECT ST_Union(ST_Buffer(ST_MakeLine(ST_LineInterpolatePoint(a.geom, 0.9), b.geom),0.00001)) geom FROM tblc a CROSS JOIN tbla b),
tble AS (SELECT ((ST_Dump(ST_DelaunayTriangles(ST_Collect(geom)))).geom) geom FROM tblb)
SELECT ST_Union(a.geom) geom FROM tble a JOIN tbld b ON ST_Intersects(a.geom, b.geom);

See the figure below for the result enter image description here

Original solutions.

The script is called: ST_PointsInStarPolygonFromDelaunayTriangles

Memo: As long as there is no universal geoinstrument that solves any geometric problems, the geo-space is much more complicated 🙂...

OR

WITH pts(geom) AS (VALUES
(st_transform(st_setsrid(st_makepoint(-97.5660461, 30.4894905), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5657216, 30.4902173), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5608779, 30.4896142), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5605001, 30.491422), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5588115, 30.4911697), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5588262, 30.4910204), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5588262, 30.4910204), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5585742, 30.4909966), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5578045, 30.4909263), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5574653, 30.4908877), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5571534, 30.4908375), 4326),4269)),
(st_transform(st_setsrid(st_makepoint(-97.5560964, 30.4907427), 4326),4269))), 
tbla AS (SELECT ST_Centroid(ST_Collect(geom)) geom FROM pts),
tblb AS (SELECT ST_MakeLine(a.geom, b.geom) geom FROM tbla a CROSS JOIN pts b),
tblc AS (SELECT ST_Union(ST_Buffer(ST_MakeLine(ST_LineInterpolatePoint(a.geom, 0.9), b.geom),0.00001)) geom FROM tblb a CROSS JOIN tbla b),
tbld AS (SELECT ((ST_Dump(ST_DelaunayTriangles(ST_Collect(geom)))).geom) geom FROM pts)
SELECT ST_Union(a.geom) geom FROM tbld a JOIN tblc b ON ST_Intersects(a.geom, b.geom);

for your geodata set, the function will look like this:

CREATE OR REPLACE FUNCTION ST_PointsInStarPolygonFromDelaunayTriangles(
geom GEOMETRY
)
RETURNS GEOMETRY AS
$BODY$
WITH
tbla AS (SELECT (ST_DumpPoints(geom)).geom geom),
tblb AS (SELECT ST_Centroid(ST_Collect(geom)) geom FROM tbla),
tblc AS (SELECT geom FROM tbla UNION SELECT geom FROM tblb),
tbld AS (SELECT ST_MakeLine(a.geom, b.geom) geom FROM tblb a CROSS JOIN tbla b),
tble AS (SELECT ST_Union(ST_Buffer(ST_MakeLine(ST_LineInterpolatePoint(a.geom, 0.9), b.geom),0.00001)) geom FROM tbld a CROSS JOIN tblb b),
tblf AS (SELECT ((ST_Dump(ST_DelaunayTriangles(ST_Collect(geom)))).geom) geom FROM tblc)
SELECT ST_Union(a.geom) geom FROM tblf a JOIN tble b ON ST_Intersects(a.geom, b.geom)
$BODY$
LANGUAGE SQL

Run:

SELECT ST_PointsInStarPolygonFromDelaunayTriangles(geom) geom FROM (SELECT ST_Collect(geom) AS geom FROM pts) a;
| improve this answer | |

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.