Missing line segment with self-intersection and MVT

[mhkeller]

I'm trying to generate a vector tile layer of a linestring consisting of four points (you could use more than four points but this is just a simplified example). Here is the linestring in geojson (​https://gist.github.com/mhkeller/9a65553c28d5063b45817b7d996b59c3)

The second and the fourth points are the same, so the line doubles back on itself. When I use this query (​https://gist.github.com/mhkeller/05d7eafc43966a2b3b5b2d41784b300b) in a tile server, the third point in the line string is missing.

I've made a full but simple reproduction with photos here: ​https://github.com/mhkeller/postgis-line-repro

I also tried using the dirt-simple-postgis-http-api (​https://github.com/tobinbradley/dirt-simple-postgis-http-api) instead of the one I wrote and I get the same result.

However, if I generate the tiles using geojson2mvt (​https://github.com/NYCPlanning/labs-geojson2mvt), the line renders correctly. A stack exchange user reported that generating the tiles with MapTiler also renders the line correctly (​https://gis.stackexchange.com/questions/471851/missing-line-segment-when-serving-a-vector-tile-from-postgis#comment772175_471851).

[pramsey]

In what XYZ tile address do you see the effect? All at all Z scales? Multiple Z scales? One? A useful repro from my PoV is a piece of SQL that generates an obviously wrong tile.

[pramsey]

For example, I tried these tests and didn't repro.

CREATE TABLE path_line AS 
SELECT 'LINESTRING(-74.00257 40.74101, -73.98698 40.74101, -73.97476 40.729633, -73.98698 40.74101)'::geometry(Linestring, 4326) AS geom;

SELECT ST_AsText(ST_AsMVTGeom(ST_Transform(geom,3857), ST_TileEnvelope(17,38598,49267))) FROM path_line;
SELECT ST_AsText(ST_AsMVTGeom(ST_Transform(geom,3857), ST_TileEnvelope(17,38597,49267))) FROM path_line;

SELECT ST_AsText(ST_AsMVTGeom(ST_Transform(geom,3857), ST_TileEnvelope(16,19298,24633))) FROM path_line;
SELECT ST_AsText(ST_AsMVTGeom(ST_Transform(geom,3857), ST_TileEnvelope(16,19299,24633))) FROM path_line;

[tjcaverly]

This issue is also described here ​https://gis.stackexchange.com/questions/471851/missing-line-segment-when-serving-a-vector-tile-from-postgis

The problem is in the algorithm used by ST_SIMPLIFY with tolerance = 0.

So

select st_astext(
    st_simplify(
    'LINESTRING (2 2,3 2,4 1,3 2, 4 4)'
    ,0)
    );

incorrectly removes the 3rd point. This will happen whenever a linestring contains consecutive points A, B and C with A=C. This will affect any vector tile geoms created with ST_ASMVTGEOM containing such linestrings, as long as the whole linestring is on one tile.

ST_SIMPLIFY tolerance = 0 uses ptarray_simplify_in_place_tolerance0, which will remove point B between A and C if it is exactly on line segment A, C, however it needs special handling for the case when A=C since this will cause false positives. This could be fixed by always keeping point B in cases with A = C. (I think removing consecutive duplicate points is handled in a separate step and doesn't need to be done here)

[robe]

when I do

select st_astext(
    st_simplify(
    'LINESTRING (2 2,3 2,4 1,3 2, 4 4)'
    ,0)
    );

I get back:

` LINESTRING(2 2,3 2,4 4)

`

Using:

` SELECT postgis_full_version();

POSTGIS="3.5.0dev 3.4.0rc1-972-ga832b8598" [EXTENSION] PGSQL="160" GEOS="3.13.0dev-CAPI-1.18.0" PROJ="9.3.1 NETWORK_ENABLED=OFF URL_ENDPOINT=​https://cdn.proj.org USER_WRITABLE_DIRECTORY=/tmp/proj DATABASE_PATH=/usr/share/proj/proj.db" GDAL="GDAL 3.8.1, released 2023/11/28" LIBXML="2.9.14" LIBJSON="0.17" LIBPROTOBUF="1.4.1" WAGYU="0.5.0 (Internal)" (core procs from "3.5.0dev 3.4.0rc1-1013-g1284f168d" need upgrade) RASTER (raster lib from "3.5.0dev 3.4.0rc1-1039-g88c3893c7" need upgrade) (raster procs from "3.5.0dev 3.4.0rc1-1013-g1284f168d" need upgrade) `

This is different from the answer stated in the stack exchange you reference of:

select st_astext(
    st_simplify(
    'LINESTRING (2 2,3 2,4 1,3 2)'
    ,0)
    );

"LINESTRING(2 2,3 2,3 2)"

What exactly are you getting? and please state the versions of postgis and related libs you are using.

[tjcaverly]

SELECT postgis_full_version()

POSTGIS="3.3.2 4975da8" [EXTENSION] PGSQL="140" GEOS="3.9.0-CAPI-1.16.2" PROJ="7.2.1" LIBXML="2.9.10" LIBJSON="0.15" LIBPROTOBUF="1.3.3" WAGYU="0.5.0 (Internal)"

I get the same result as you,

select st_astext(
    st_simplify(
    'LINESTRING (2 2,3 2,4 1,3 2, 4 4)'
    ,0)
    );

"LINESTRING(2 2,3 2,4 4)"

This incorrectly removes the 3rd point, erasing a large portion of the geometry. It then removes the 4th point since it's the same as the 2nd, which doesn't affect anything. Correct behavior would be returning LINESTRING (2 2,3 2,4 1,3 2, 4 4).

The example on stackexchange is different but shows the same bug.

[pramsey]

Hm, manifests in a fun performance optimization, and gives the right answer if that optimization is removed.

​https://github.com/postgis/postgis/blame/1d12c312bdb2b7d900a01d0259b1f1259edf211d/liblwgeom/ptarray.c#L1667

[pramsey]

This fixed the particular case, and passes all tests ​https://github.com/pramsey/postgis/commit/8c30a25acb166de09a7d647e1e02dc6d71eae9ea

diff --git a/liblwgeom/ptarray.c b/liblwgeom/ptarray.c
index 7d9c8c0cd..b331c9558 100644
--- a/liblwgeom/ptarray.c
+++ b/liblwgeom/ptarray.c
@@ -1685,7 +1685,11 @@ ptarray_simplify_in_place_tolerance0(POINTARRAY *pa)
                double dot_ac_ab = ca_x * ba_x + ca_y * ba_y;
                double s_numerator = ca_x * ba_y - ca_y * ba_x;
 
-               if (dot_ac_ab < 0.0 || dot_ac_ab > ab_length_sqr || s_numerator != 0)
+               if (p2d_same(kept_pt, next_pt) ||
+                       dot_ac_ab < 0.0 ||
+                       dot_ac_ab > ab_length_sqr ||
+                       s_numerator != 0)
+
                {
                        kept_it++;
                        kept_pt = curr_pt;

[Paul Ramsey <pramsey@…>]

In 2c68d130/git:

Do not remove points when doubling back at zero tolerance, references #5654

[Paul Ramsey <pramsey@…>]

In 8c30a25a/git:

Do not remove points when doubling back at zero tolerance, references #5654

[pramsey]

Fixed at [8c30a25ac/git]