Closing a concave V-shaped MultiLineString in PostGIS to form a Polygon

Question

I have a 3D MultiLineString (roughly ~2'400 segments) in PostGIS forming kind of a V-shape (blue line on the image bellow) which I need to close (e.g. with the desired green closing straight line) in order to form a Polygon:

I found ST_MakePolygon which looked promising, but it raises an error:

ERROR:  Shell is not a line
SQL state: XX000

This is certainly because the line is actually of type ST_MultiLineString.

So I tried to figure out a way to convert this to a true LineString but the usage of ST_LineMerge doesn't seem to have any effect on the resulting geometry type:

SELECT ST_GeometryType(ST_LineMerge(geom)) AS geomtype;
  geomtype      


ST_MultiLineString
(1 row)

Any ideas on how to succeed ?
The line should be connected, but based on these observations, it's apparently not.
The differences (i.e. gaps?), if existing, should be due to numerical approximations.
The original geometry is valid (ST_IsValid returns true).

EDIT:

I confirm the gaps (in the blue multilinestring) due to computational issues, e.g. with these two lines which should be connected with their first respective points:

MULTILINESTRING Z (( 2702603.3370974000 1253383.3314513000 607.1897036409243, <POINT_B> ))
MULTILINESTRING Z (( 2702603.3370974010 1253383.3314512996 607.1897036407847, <POINT_D> ))

Those are two original parts which are aggregated to actually compose the blue multilinestring.

Note: I appended trailing zeroes in coordinates to keep numbers vertically aligned for a better visual perception but they were not present in the query results.

EDIT 2:

ST_MakePolygon() doesn't seem to work on a mathematically perfectly connected multilinestring, but it actually throws another error:

SELECT
ST_MakePolygon(
    ST_LineMerge(
        ST_GeomFromText(
            'MULTILINESTRING(
            (2702703.3370974000 1253383.3314513000 607.1897036409243,
             2702779.0888727466 1253368.1777454934 596.9877159267046),
            (2702703.3370974000 1253383.3314513000 607.1897036409243,
             2702638.0847548526 1253396.2982941763 613.5631635554059),
            (2702638.0847548526 1253396.2982941763 613.5631635554059,
             2702678.5554887933 1253321.3220054788 609.5533124641558),
            (2702678.5554887933 1253321.3220054788 609.5533124641558,
             2702682.5996487933 1253250.0846754788 601.2201976641112)
            )',
            2056
        )
    )
)

This raises:

ERROR:  lwpoly_from_lwlines: shell must be closed
SQL state: XX000

And the ST_ConcaveHull(), as suggested by @hbg in the comment below, is not working either:

SELECT
ST_ConcaveHull(
    ST_LineMerge(
        ST_GeomFromText(
            'MULTILINESTRING(
            (2702703.3370974000 1253383.3314513000 607.1897036409243,
             2702779.0888727466 1253368.1777454934 596.9877159267046),
            (2702703.3370974000 1253383.3314513000 607.1897036409243,
             2702638.0847548526 1253396.2982941763 613.5631635554059),
            (2702638.0847548526 1253396.2982941763 613.5631635554059,
             2702678.5554887933 1253321.3220054788 609.5533124641558),
            (2702678.5554887933 1253321.3220054788 609.5533124641558,
             2702682.5996487933 1253250.0846754788 601.2201976641112)
            )',
            2056
        )
    )
,0.1
)

which raises:

ERROR:  Cannot ST_Collect geometries with differing dimensionality.
CONTEXT:  PL/pgSQL function _st_concavehull(geometry) line 67 at assignment
PL/pgSQL function st_concavehull(geometry,double precision,boolean) line 152 at assignment
PL/pgSQL assignment "var_tempgeom := public.ST_Buffer(public.ST_ConcaveHull(var_geoms[1],least(param_pctconvex + param_pctconvex/var_div),true),var_buf, 'quad_segs=2')"
PL/pgSQL function st_concavehull(geometry,double precision,boolean) line 114 at assignment
PL/pgSQL assignment "var_tempgeom := public.ST_Buffer(public.ST_ConcaveHull(var_geoms[1],least(param_pctconvex + param_pctconvex/var_div),true),var_buf, 'quad_segs=2')"
PL/pgSQL function st_concavehull(geometry,double precision,boolean) line 114 at assignment
PL/pgSQL assignment "var_tempgeom := public.ST_Buffer(public.ST_ConcaveHull(var_geoms[1],least(param_pctconvex + param_pctconvex/var_div),true),var_buf, 'quad_segs=2')"
...the two previous lines being repeated multiple times.
PL/pgSQL function st_concavehull(geometry,double precision,boolean) line 114 at assignment
SQL state: XX000

Notice: the same errors on the two last examples also occur with PG 14 + PostGIS 3.3.0 (docker image: postgis/postgis:14-3.3.0rc2-alpine based on GEOS="3.10.2-CAPI-1.16.0")

PostGIS 3.3.1 with the very latest PG 15 (beta4) is still based on GEOS 3.10.2:
POSTGIS="3.3.1 0" [EXTENSION] PGSQL="150" GEOS="3.10.2-CAPI-1.16.0".

Version info:

PostgreSQL 14.0 (Debian 14.0-1.pgdg110+1) on x86_64-pc-linux-gnu, compiled by gcc (Debian 10.2.1-6) 10.2.1 20210110, 64-bit

POSTGIS="3.1.4 ded6c34" [EXTENSION] PGSQL="140" GEOS="3.9.0-CAPI-1.16.2" SFCGAL="1.3.8" PROJ="7.2.1" LIBXML="2.9.10" LIBJSON="0.15" LIBPROTOBUF="1.3.3" WAGYU="0.5.0 (Internal)" TOPOLOGY

It might be nice if ST_LineMerge accepted a distanceTolerance parameter so it could snap close endpoints together. — dr_jts, Sep 09 '22 at 19:36
Indeed, especially for such micro or nanometric differences ! :'( — swiss_knight, Sep 10 '22 at 11:40
I would try concave hull https://postgis.net/docs/ST_ConcaveHull.html. — user30184, Sep 10 '22 at 13:18
You've shown two multilinestrings (3D) above. Do those include the green line in your diagram or is that separate? The dimensionality error suggests to me that you are mixing 3D (multilinestring Z) with a 2D linestring when calling ST_ConcaveHull. Have you tried an example with identical endpoints in your multilinestring to see if ST_LineMerge and ST_MakePolygon work in that case without worrying about the snapping tolerance? — hgb, Sep 10 '22 at 17:04
The two 3D multilinestrings examples shown above are two original segments composing the blue multilinestring. The green one is the desired one. It doesn't exist yet. I will try to build an artificial data set with strictly identical points to show the desired result then, but I guess it will work. Hopefully. — swiss_knight, Sep 10 '22 at 17:14
ST_ConcaveHull has been improved in PostGIS 3.3/GEOS 3.11, so you could try that. While it may seem odd, theoretically it could provide a polygon that respects the input line if the tolerance value is small. — dr_jts, Sep 10 '22 at 18:08
If you don't need full fidelity with the original points, you could use ST_SnapToGrid to remove the excess precision causing the discrepancy between endpoints. — dr_jts, Sep 10 '22 at 18:10
I added a reproducible example with a mathematically connected multiline (exact same coordinates on the points of junction), which doesn't give me the desired results unfortunately, but with another error. I also added the ST_ConcaveHull() trick and its results. @dr_jts I also tried PostGIS 3.3.0 but the errors persist. See last update in main post. — swiss_knight, Sep 10 '22 at 18:28
ST_MakePolygon does require rings to be closed (have duplicate start and end points. — dr_jts, Sep 10 '22 at 19:04
The PG 3.3 Docker image is using GEOS 3.10, so still uses the old concave hull algorithm. — dr_jts, Sep 10 '22 at 19:04

score 2 · Answer 1 · answered Sep 11 '22 at 15:32

Building on your ST_MakePolygon example. As @dr_jts pointed out in a comment, ST_MakePolygon won't work unless you give it a closed ring. Here's a crude consrtuctor, relying on your statement that all input geometries form a "V".

WITH merge AS (
  SELECT ST_LineMerge(
    ST_GeomFromText(
      'MULTILINESTRING(
      (2702703.3370974000 1253383.3314513000 607.1897036409243,
      2702779.0888727466 1253368.1777454934 596.9877159267046),
      (2702703.3370974000 1253383.3314513000 607.1897036409243,
      2702638.0847548526 1253396.2982941763 613.5631635554059),
      (2702638.0847548526 1253396.2982941763 613.5631635554059,
      2702678.5554887933 1253321.3220054788 609.5533124641558),
      (2702678.5554887933 1253321.3220054788 609.5533124641558,
      2702682.5996487933 1253250.0846754788 601.2201976641112)
      )',
      2056
    )
  ) AS line
)
SELECT ST_MakePolygon(ST_MakeLine(line, ST_PointN(line, 1))) AS geom
  FROM merge;

Line merge from the multilinestring, as you showed in your example
Explicitly copy the first point to the end of the linestring to create a closed loop. Depending on your geometries, this could create an invalid geometry.
Make polygon from the created ring.

In your real geometries, you still need to deal with the gaps, possibly snapping to grid as also suggested by @dr_jts.

I accept this because it was posted first but @dr_jts' answer is more especially suited in the case of numerical inaccuracies. — swiss_knight, Oct 02 '22 at 14:41

score 1 · Answer 2 · answered Sep 14 '22 at 00:24

1

If the line elements are in the MultiLineString, you can use ST_Dump and ST_MakeLine to "sew" them together into a single line, as in this answer: https://gis.stackexchange.com/a/375977/14766. Then append a closing point and use ST_MakePolygon.

answered Sep 14 '22 at 00:24

dr_jts

5,038
18
15

Your answer is best suited for solving the issue in case of numerical inaccuracies. I've keep things separated in a new thread. See: https://gis.stackexchange.com/q/441991/65370 where I may benefit from some simplifications, I guess. – swiss_knight Oct 02 '22 at 14:41

score 0 · Answer 3 · answered Oct 05 '22 at 18:47

I'm not sure I fully understand your question, but with the geodata sample provided, and the proposed Geo pgSQL architecture, everything is fine for me.

So one more option.

Input:

WITH 
pts(geom) AS (VALUES (ST_GeomFromText('MULTILINESTRING(
            (2702703.3370974000 1253383.3314513000 607.1897036409243,
             2702779.0888727466 1253368.1777454934 596.9877159267046),
            (2702703.3370974000 1253383.3314513000 607.1897036409243,
             2702638.0847548526 1253396.2982941763 613.5631635554059),
            (2702638.0847548526 1253396.2982941763 613.5631635554059,
             2702678.5554887933 1253321.3220054788 609.5533124641558),
            (2702678.5554887933 1253321.3220054788 609.5533124641558,
             2702682.5996487933 1253250.0846754788 601.2201976641112)
            )', 2056))),
tbla AS (SELECT ST_Force2D((ST_DumpPoints(geom)).geom) geom FROM pts),
centroid AS (SELECT ST_Centroid(ST_Collect(geom)) centroid FROM pts),
tblb AS (SELECT ST_MakeLine(geom ORDER BY ST_Azimuth(centroid, geom)) geom FROM tbla CROSS JOIN centroid)
SELECT ST_MakePolygon(ST_MakeLine(geom, ST_PointN(geom, 1))) AS geom FROM tblb

Output:

Picture

FOGS

SELECT ST_MakePolygon(ST_AddPoint(geom, ST_StartPoint(geom))) geom FROM tblb — Cyril Mikhalchenko, Oct 05 '22 at 19:24

Closing a concave V-shaped MultiLineString in PostGIS to form a Polygon

EDIT:

EDIT 2:

3 Answers3

Linked