Bug in geography ST_Segmentize

[robe]

There is still something wrong with ST_Segmentize geography as someone noted on IRC:

From ​http://irclogs.geoapt.com/postgis/%23postgis.2016-11-08.log

15:42:57	skyrocker:	E.g. SELECT ST_AsText(ST_Segmentize(ST_GeographyFromText('LINESTRING(38.769917 10.780694, 38.769917 9.106194)'), 30000))
15:43:01	skyrocker:	Results in LINESTRING(38.769917 10.780694,3.97245555092463e-308 5.68175492717434e-322,38.769917 9.106194)
15:43:12	skyrocker:	Which is obviously not correct
15:45:18	skyrocker:	Although after changing first coordinate of first point in LINESTRING I get correct result (e.g. instead of 38.769917 put 38.769918)
15:46:50	skyrocker:	And also for some werid reason swapping order of points in LINESTRING also gives me correct result

The only thing I see special about this is that in cartesian it's vertical.

Testing on

POSTGIS="2.3.0 r15146" GEOS="3.6.0-CAPI-1.10.0 r4265" PROJ="Rel. 4.9.1, 04 March 2015" GDAL="GDAL 2.1.1, released 2016/07/07" LIBXML="2.7.8" LIBJSON="0.12" RASTER

The lengths are way off, so it's clearly wrong:

SELECT ST_Length(geog) As len_before, ST_Length(ST_Segmentize(geog, 30000)) As len_after
FROM (SELECT ST_GeographyFromText('LINESTRING(38.769917 10.780694, 38.769917 9.106194)') As geog) AS f;

    len_before    |    len_after
------------------+------------------
 185212.200983372 | 8865756.13359671
(1 row)

But as noted flipping the coordinates works

SELECT ST_Length(geog) As len_before, ST_Length(ST_Segmentize(geog, 30000)) As len_after
FROM (SELECT ST_GeographyFromText('LINESTRING(38.769917 9.106194, 38.769917 10.780694)') As geog) AS f;

    len_before    |    len_after
------------------+------------------
 185212.200983372 | 185212.200983372
(1 row)

This problem is new in 2.3.0 but not 2.2.2, so must be a continuation of issue I complained about in #3539:

In POSTGIS="2.2.2 r14797" GEOS="3.5.0-CAPI-1.9.0 r4090" PROJ="Rel. 4.9.1, 04 March 2015" GDAL="GDAL 2.0.2, released 2016/01/26" LIBXML="2.7.8" LIBJSON="0.12" RASTER the answers are right:

SELECT ST_Length(geog) As len_before, ST_Length(ST_Segmentize(geog, 30000)) As len_after
FROM (SELECT ST_GeographyFromText('LINESTRING(38.769917 10.780694, 38.769917 9.106194)') As geog) AS f;

    len_before    |    len_after
------------------+------------------
 185212.200983372 | 185212.200983372
(1 row)

[hmercier]

Should be fixed with ​https://github.com/postgis/postgis/pull/119

[robe]

In 15222:

Bug in geography ST_Segmentize
References #3667 for PostGIS 2.3.1
Patch provided by Hugo Mercier (Oslandia)

[robe]

In 15223:

Bug in geography ST_Segmentize
Closes #3667 for PostGIS 2.4.0 (trunk)
Patch provided by Hugo Mercier (Oslandia)

[robe]

Someone reported on IRC that this is still not fixed with this example. I have confirmed it is still broken.

I know it's a breaking change, but can we revert back to the old pre 2.3 behavior at least for PostGIS 2.4.1 (2.3.4 would be nice to do too, though I'm okay with letting that one go).

SELECT ST_AsText(ST_Segmentize(ST_GeogFromText('LINESTRING(11.3953811 47.2577728,11.3953811 47.257936)'), 15));

outputs

                                   st_astext
--------------------------------------------------------------------------------
 LINESTRING(11.3953811 47.2577728,0 5.68175492717434e-322,11.3953811 47.257936)
(1 row)

Which is huge. The orginal line was only 19 meters long. I also tested on 2.3.3 and got the same stupid answer.

The old behavior may not have been as good about splitting the lines evenly, but at least it always gave reasonable answers.

IRC thread from Nikolaus Krismer

13:25:25	Niko_K:	Can I please get some help on the ST_Segmentize function in Postgis 2.4.0?
13:25:47	Niko_K:	I think that the bug #3667 (https://trac.osgeo.org/postgis/ticket/3667) is not fixed
13:25:48	sigq:	Title: #3667 (Bug in geography ST_Segmentize) – PostGIS (at trac.osgeo.org)
13:28:14	Niko_K:	Just try the following SQL command
13:28:16	Niko_K:	SELECT ST_AsText(ST_Segmentize(ST_GeogFromText('LINESTRING(11.3953811 47.2577728,11.3953811 47.257936)'), 15))
13:28:26	Niko_K:	For me the result is:
13:28:28	Niko_K:	LINESTRING(11.3953811 47.2577728,8.58553798490362e-315 14376.1464672296,11.3953811 47.257936)
13:29:11	Niko_K:	the inserted point '8.58553798490362e-315 14376.1464672296' is obviously not correct
13:29:17	Niko_K:	some information on the version used
13:29:25	Niko_K:	POSTGIS="2.4.0 r15853" PGSQL="100" GEOS="3.6.1-CAPI-1.10.1 r0" PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 2.1.3, released 2017/20/01" LIBXML="2.9.4" LIBJSON="0.12.1"
13:29:27	Niko_K:	LIBPROTOBUF="1.2.1" RASTER
13:30:07	Niko_K:	System is fedora 26 with Postgis installed from latest PostgreSQL repo
13:30:53	Niko_K:	PostgreSQL version is: "PostgreSQL 10.0 on x86_64-pc-linux-gnu, compiled by gcc (GCC) 7.2.1 20170915 (Red Hat 7.2.1-2), 64-bit"
13:45:34	Niko_K:	Just found out that the ST_Segmentize sometimes leads to intermediate points, with lat and lon close to zero and some (like the example above) where only one of lat or lon is close to zero
13:45:45	Niko_K:	here lat and lon become close to zero:
13:45:47	Niko_K:	SELECT ST_AsText(ST_Segmentize(ST_GeomFromText('LINESTRING(11.4204246 47.2494256,11.4203095 47.2493851,11.4203095 47.2492321,11.4205452 47.2489215)', 4326)::geography, 15))
13:46:17	Niko_K:	in the example above only one is close to zero while to other one is too large (larger than 180 and 90 degrees)
13:46:22	Niko_K:	SELECT ST_AsText(ST_Segmentize(ST_GeogFromText('LINESTRING(11.3953811 47.2577728,11.3953811 47.257936)'), 15))

[nikolauskrismer]

User creation took some time… sorry to be a bit late here
If you need more information about my setup or need more test cases, I am happy to help

[pramsey]

I have an approach which hopefully will retain the simplicity of the original while returning a segmentization made of all equal lengths.

[pramsey]

In 16092:

Replace project-and-extend with bisect-and-recurse method
for generating segmentized geography. Provides "mostly equal"
segment lengths but also has more numerical stability
for small cases, as the old 3-space projection approach
did. Closes #3667

[pramsey]

I closed this with a fix to trunk. But, it should perhaps come back to 2.3/2.4. That would involve a behaviour change though, since the exact segmentizations returned are not the same. Thoughts?

[robe]

Just do it. I don't think most people noticed the behavior improvement in the first place and the introduced bug corner cases are definitely noticeable and annoying.

[pramsey]

In 16094:

Replace project-and-entend logic with
bisect-and-recurse in geography segmentization.
Preserves "mostly equal" segment lengths, and
should be more numerically stable.
Backport to 2.4.
References #3667

[pramsey]

In 16095:

Replace project-and-entend logic with
bisect-and-recurse in geography segmentization.
Preserves "mostly equal" segment lengths, and
should be more numerically stable.
Backport to 2.3.
References #3667

[robe]

In 16098:

Add missing release notes for ST_Segmentize geography fix.
References #3667 for PostGIS 2.4.2

[robe]

In 16100:

Add missing release notes for ST_Segmentize geography fix.
References #3667 for PostGIS 2.3.5