Opened 4 years ago

Closed 3 years ago

Last modified 3 years ago

#3667 closed defect (fixed)

Bug in geography ST_Segmentize

Reported by: robe Owned by: pramsey
Priority: blocker Milestone: PostGIS 2.3.5
Component: postgis Version: 2.3.x
Keywords: Cc: hmercier

Description

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)

Change History (15)

comment:1 Changed 4 years ago by robe

Cc: hmercier added

comment:2 Changed 4 years ago by hmercier

comment:3 Changed 4 years ago by robe

In 15222:

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

comment:4 Changed 4 years ago by robe

Resolution: fixed
Status: newclosed

In 15223:

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

comment:5 Changed 3 years ago by robe

Milestone: PostGIS 2.3.1PostGIS 2.3.4
Priority: highblocker
Resolution: fixed
Status: closedreopened

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))
Last edited 3 years ago by robe (previous) (diff)

comment:6 Changed 3 years ago by 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

comment:7 Changed 3 years ago by pramsey

Milestone: PostGIS 2.3.4PostGIS 2.3.5

comment:8 Changed 3 years ago by pramsey

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

comment:9 Changed 3 years ago by pramsey

Resolution: fixed
Status: reopenedclosed

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

comment:10 Changed 3 years ago by 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?

comment:11 Changed 3 years ago by 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.

comment:12 Changed 3 years ago by 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

comment:13 Changed 3 years ago by 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

comment:14 Changed 3 years ago by robe

In 16098:

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

comment:15 Changed 3 years ago by robe

In 16100:

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

Note: See TracTickets for help on using tickets.