CODEWITHSUNDEEP
postgispostgresql-9.3
djkern
djkern

Reputation: 482

Unexplained behavior of st_linesubstring

Given the documentation of ST_LINESUBSTRING (https://postgis.net/docs/ST_LineSubstring.html), you would expect that it would be straightforward to divide a 2d linestring into two equal parts.

See below:

SELECT ST_LENGTH( ST_LINESUBSTRING( geom, 0.0, 0.5 )::GEOGRAPHY ),
       ST_LENGTH( ST_LINESUBSTRING( geom, 0.5, 1.0 )::GEOGRAPHY )
FROM data.street WHERE id = '11242';

    st_length     |    st_length
------------------+------------------
2003.80912128768 | 2153.56113185375 (1 row)

Notice that the command executed above should have yielded two halves of equal length. However, the halves are of differing lengths.

Given the following, you would expect each half to be ~2078:

SELECT ST_LENGTH( geom::GEOGRAPHY )
FROM data.street WHERE id = '11242';
    st_length    
-----------------
4157.37025314031
(1 row)

At least, if you add each half together, you get 4157.370253, which is the correct length. Thus, ST_LINESUBSTRING isn't shrinking or expanding the geometry.

So, does anyone have an answer? Why does ST_LINESUBSTRING behave this way?

Upvotes: 2

Views: 689

Answers (1)

ewcz
ewcz

Reputation: 13097

The reason for this behavior is that ST_LineSubstring splits geometries, not geographies. In other words, it does the splitting in the projection coordinates in which the LineStrings are stored.

For example, let's assume that a particular line connects two points A, B on the globe. Now, if one projects these points into lon/lat coordinates and calculates the midpoint, the result will be different than projecting the midpoint of the great circle connecting A and B.

One solution would be to use ST_Segmentize:

CREATE TEMPORARY TABLE lines(
  geom GEOMETRY
);

INSERT INTO lines VALUES
  (ST_MakeLine(ST_MakePoint(-13.4, 52.25), ST_MakePoint(-11.582, 48.1351)))
;

SELECT ST_Length(geom::geography) FROM lines;

SELECT
  ST_Length( ST_LineSubstring(geom, 0.0, 0.5 )::GEOGRAPHY ),
  ST_Length( ST_LineSubstring(geom, 0.5, 1.0 )::GEOGRAPHY )
FROM lines;

SELECT
  ST_AsText(ST_MakeLine(start_point, end_point)) AS segment,
  ST_Length(ST_MakeLine(start_point, end_point)::geography) AS len
FROM
(
    SELECT
        ST_Pointn(geom, generate_series(1, ST_NumPoints(geom)-1)) as start_point,
        ST_Pointn(geom, generate_series(2, ST_NumPoints(geom))) as end_point
    FROM (
        SELECT ST_Segmentize(geom::geography, ST_Length(geom::geography)/2)::geometry AS geom FROM lines
    ) AS line
) AS tmp;

this produces:

CREATE TABLE
INSERT 0 1
    st_length     
------------------
 475724.617118868
(1 row)

    st_length    |    st_length     
-----------------+------------------
 237537.46004128 | 238220.186078105
(1 row)

                            segment                             |       len        
----------------------------------------------------------------+------------------
 LINESTRING(-13.4 52.25,-12.4518126471826 50.1960897877664)     |  237904.63855338
 LINESTRING(-12.4518126471826 50.1960897877664,-11.582 48.1351) | 237819.978632575
(2 rows)

Using the latter approach, the segments are indeed of approximately equal length.

Upvotes: 4

Related Questions