ST_LineSubstring returns incorrect lengths in geography

[jorsu]

I am encountering an issue with the ST_LineSubstring function. The results remain inconsistent when applying on a geographytype.

I am running:

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

POSTGIS="3.5.2 dea6d0a" [EXTENSION] PGSQL="170" GEOS="3.9.0-CAPI-1.16.2" PROJ="7.2.1 NETWORK_ENABLED=OFF URL_ENDPOINT=https://cdn.proj.org USER_WRITABLE_DIRECTORY=/var/lib/postgresql/.local/share/proj DATABASE_PATH=/usr/share/proj/proj.db" (compiled against PROJ 7.2.1) GDAL="GDAL 3.2.2, released 2021/03/05" LIBXML="2.9.10" LIBJSON="0.15" LIBPROTOBUF="1.3.3" WAGYU="0.5.0 (Internal)" TOPOLOGY RASTER

I am trying to:

• Extract a subline from a LINESTRING.
• Accurately measure its length in the geography type.

However, when using ST_LineSubstring directly on a geography type, the resulting length is incorrect and does not match the expected value.

Here is a minimal reproducible example:

WITH
    "path" AS (
        SELECT ST_SetSRID(ST_GeomFromText('LINESTRING(2.3522 48.8566, 2.3535 48.8570, 2.3548 48.8580, 2.3561 48.8590, 2.3574 48.8600)'), 4326) AS path
    )
SELECT
    ST_Length(p."path"::geography) full_length,
    ST_Length(
         ST_LineSubstring(p."path"::geometry, 0, 1)::geography
    ) geom_length,
    ST_Length(
         ST_LineSubstring(p."path"::geography, 0, 1)::geography
    ) geog_length
from "path" p;

This query returns:

   full_length 	       geom_length 	   geog_length
544.8029520664281   544.8029520664281 	324.1448393269909

Why does ST_LineSubstring on a geography type produce a different length than the same operation on a geometry type and the original path?

[jratike80]

Originally asked on gis.stackexchange ​https://gis.stackexchange.com/questions/491461/st-linesubstring-returns-incorrect-lengths-in-geography-in-postgis/491468#491468.

It seems that the problem can be simplified into question "why ST_LineSubstring with range 0,1 in the following case drops the last vertex?".

SELECT
ST_Astext(
ST_LineSubstring(ST_SetSRID(ST_GeomFromText('LINESTRING(2.3522 48.8566, 2.3535 48.8570, 2.3548 48.8580, 2.3561 48.8590, 2.3574 48.8600)'), 4326)
::geography,0,1))

"LINESTRING(2.3522 48.8566,2.3535 48.857,2.3548 48.858,2.355442128456938 48.85849395179974)"

[pramsey]

Simplified version of this is even more crazy-making.

SELECT
 ST_Astext(
  ST_LineSubstring(
   'LINESTRING(2 48, 3 49, 4 50)'::geography, 0, 1));

LINESTRING(2 48,3 49,3.493455272486287 49.49965167325753)

[Paul Ramsey <pramsey@…>]

In 0f5ece1/git:

Don't skip the last point in doing
linestring substring for geography,
closes #5874

[Paul Ramsey <pramsey@…>]

In b5f353e/git:

Don't skip the last point in doing
linestring substring for geography,
References #5874

[Paul Ramsey <pramsey@…>]

In a10ca1a9/git:

Don't skip the last point in doing
linestring substring for geography,
references #5874