Bug in geography ST_Segmentize

[petermj]

I know #3667 has fixed a bug in ST_Segmentize, but I am still seeing issues.

I am segmenting a line, and depending upon the order of the points it changes the number of segments produced, with the final segment being 'massive'.

I'm sorry if my SQL code is horrendous, I am just learning.

I have/will attached an SQL query that highlights the bug.

[petermj]

SQL query test case.

[petermj]

Sorry, I should add I am running PostgreSQL 9.6.2-2, and PostGIS 2.3.2, both x64 bit on Windows 10. I have downloaded the EDB compiled version if that helps?.

[robe]

Okay I see your point. I reduced the query to one offending example:

IN POSTGIS="2.3.2 r15302" GEOS="3.6.1-CAPI-1.10.1 r4317" PROJ="Rel. 4.9.1, 04 March 2015" GDAL="GDAL 2.1.3, released 2017/20/01" LIBXML="2.7.8" LIBJSON="0.12" RASTER

With f AS  (SELECT ST_Segmentize(ST_GeogFromText('SRID=4326; LINESTRING(-4.0 49.8, -4.0 59.0)'),1000) AS geog_s, ST_GeogFromText('SRID=4326; LINESTRING(-4.0 49.8, -4.0 59.0)') AS geog),
    f2 AS (SELECT ST_NPoints(geog_s::geometry) AS geog_s_npoints, geog, geog_s FROM F)
SELECT 	geog_s_npoints, ST_Length(geog_s) AS seg_length, ST_Length(geog) AS orig_length, ST_Length(geog)/1000 As exp_np,
        ST_Length(ST_MakeLine(ST_PointN(geog_s::geometry, 1), ST_PointN(geog_s::geometry, 2))::geography) As first_seg_length,
        ST_Length(ST_MakeLine(ST_PointN(geog_s::geometry, geog_s_npoints - 1), ST_PointN(geog_s::geometry, geog_s_npoints))::geography) As last_seg_length
FROM f2;

Yields:

 geog_s_npoints |    seg_length    |   orig_length   |     exp_np      | first_seg_length | last_seg_length
----------------+------------------+-----------------+-----------------+------------------+------------------
            894 | 1024067.47066481 | 1024067.4706648 | 1024.0674706648 | 1000.26657219012 | 131222.895347768
(1 row)

Answer in POSTGIS="2.2.2 r14797" GEOS="3.6.1-CAPI-1.10.1 r4317" PROJ="Rel. 4.9.1, 04 March 2015" GDAL="GDAL 2.1.3, released 2017/20/01" LIBXML="2.7.8" LIBJSON="0.12" RASTER

 geog_s_npoints |   seg_length    |   orig_length   |     exp_np      | first_seg_length | last_seg_length
----------------+-----------------+-----------------+-----------------+------------------+------------------
           1024 | 1024067.4706648 | 1024067.4706648 | 1024.0674706648 | 995.986340600015 | 997.506446092082
(1 row)

Much more pleasant indeed.

SO yes this looks like another regression resulting from the 2.3 segmentize change.

[petermj]

b.t.w. thank you for your work on this

[robe]

this is still a problem in 2.4.0dev r15774

I'm really tempted at this point to revert this whole segmentize change we did if we can't fix this. I feel the old behavior was much better.

[robe]

Since we are seriously thinking of going back to the old behavior which was a cleaner solution (less failure modes). Have to do this in 2.4.0

[robe]

Note brought up on IRC and Mike Toews responded.

https://lists.osgeo.org/pipermail/postgis-devel/2017-September/026475.html

On 22 September 2017 at 07:47, Regina Obe <lr at pcorp.us> wrote:
> Except for this last one which I propose we fix by switching back to 2.2
> behavior.  The new behavior is harder to get right in all cases and I don't
> see it as being that much of an improvement over the old behavior.

I'll admit that I haven't been following the old or new behaviours of
ST_Segmentize, but I'd suggest redesigning this based on GeographicLib
(same library used for ST_Distance, via PROJ), since it has reliable
methods to do parts of the technique. In particular, it will produce
equal-length sub-segments, probably within a few nanometres.

Charles has a few examples of interpolating points from an inverse
geodesic line, here is one for C:
https://geographiclib.sourceforge.io/html/C/geodesic_8h.html#a0b91285fb3cd572d2d77b2e062e8c800

Essentially, you create an geod_geodesicline for each vertex pair,
find the geodesic length (s12). If the geodesic segment length is
greater than the segmentize distance, then interpolate the points
along the geodesic line with geod_position by dividing s12 by the
appropriate number of points required. Each sub-segment would have
equal distance between them.

As for performance, creating each geod_geodesicline is just as
expensive as ST_Distance (2.62 µs), then geod_position to get the
interpolated points is very cheap (each 0.37 µs).

Cheers, Mike

with that thought in mind I guess we could live with this bug for 2.4 and in 2.5 reimplement with GeographicLib.

Though that begs the question should we then in 2.5 require proj 4.9 as a minimum (or we put guards in place for lower proj).

[robe]

if we are going to switch behavior back in 2.5, might not be worth fixing in 2.4. better off with a rewrite in 2.5 and proj 4.9 as minimum.

[pramsey]

How about the old (pre-2.3) behaviour if we have proj < 4.9 and new fancy behaviour if proj >= 4.9.

[robe]

I'm good with that.

[pramsey]

This has been superceded by the fix in #3667