ST_Intersection broken with Z coordinates

[TBL]

ST_Intersection doesn't handle Z coordinates but doesn't fail when trying to do so. Instead it just returns broken data.

This shoud either be documented or (better) ST_Intersection should always drop the Z or (best) ST_Intersection should fail with a error.

[pramsey]

You clearly have a platonic idea of what ST Intersection should do in your head, but it does indeed "handle" Z, it just does it in a particular way. Basically as 2.5D coordinates. So the intersection always happens in the X/Y plane, and any new introduced vertices are given averaged values based on the inputs. It's ugly, but it was judged to be "better than nothing", which is different from your viewpoint ("do nothing!").

Unless, you have a real example of wrongness, in which case, please show some examples that demonstrate it.

[TBL]

I just posted that on freenode:

with linestring as (select ST_GeomFromText('LINESTRING Z (-1 -1 0, -0.5 -0.5 1, 0 0 2, 0.5 0.5 3, 1 1 4, 1.5 1.5 5, 2 2 6)', 4326) as geom),
     polygon as (select ST_GeomFromText('POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))', 4326) as geom)
select st_astext(st_intersection(linestring.geom, polygon.geom)) from linestring, polygon

Result is as expected: "LINESTRING Z (0 0 2,0.5 0.5 3,1 1 4)"

But with a slightly longer linestring:

with linestring as (select ST_GeomFromText('LINESTRING Z (-20.00 -20.00 0, -19.95 -19.95 1, [...], 19.990 19.990 7998, 19.995 19.995 7999)', 4326) as geom),
     polygon as (select ST_GeomFromText('POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))', 4326) as geom)
select st_astext(st_intersection(linestring.geom, polygon.geom)) from linestring, polygon

"MULTILINESTRING Z ((0 0 400,0.005 0.005 2200.55),(0.005 0.005 2200.55,0.01 0.01 2201.1),(0.01 0.01 2201.1,0.015 0.015 2201.65),(0.015 0.015 2201.65,0.02 0.02 2202.2),(0.02 0.02 2202.2,0.025 0.025 2202.75),(0.025 0.025 2202.75,0.03 0.03 2203.3),(0.03 0.03 2 (…)"

The z coordinates of the result are somehow mixed up, while the x/y coordinates seem to be right (the z at 0 0 should be 4000, at 0.005 0.005 it should be 4001 …). The returned values don't make any sense in my eyes.

[pramsey]

Since I don't have your exact input line, I cannot duplicate, but building a regular line programmatically, I'm not seeing the same result:

with linestring as (
select st_setsrid(st_makeline(st_makepoint((a-5000)/1000.0,(a-5000)/1000.0,a)),4326) as geom from generate_series(0,10000) as a
),
polygon as (
select ST_GeomFromText('POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))', 4326) as geom)
select st_astext(st_intersection(linestring.geom, polygon.geom)) from linestring, polygon

 POSTGIS="2.1.4dev r12671" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0, 6 March 2012" LIBXML="2.9.0" LIBJSON="UNKNOWN"

[TBL]

Sorry, I fucked up the example. Try this one:

with linestring as (
  select st_makeline(
    (select st_setsrid(st_makeline(st_makepoint((a-400)/20.0,(a-400)/20.0,a)),4326) from generate_series(0,800) as a),
    (select st_setsrid(st_makeline(st_makepoint(-(a-1200)/20.0,-(a-1200)/20.0,a)),4326) from generate_series(801,1600) as a)
  ) as geom
),
polygon as (select ST_GeomFromText('POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))', 4326) as geom)
select st_astext(st_intersection(linestring.geom, polygon.geom)) from linestring, polygon

Expected output would be

MULTILINESTRING Z ((0 0 400, 0.05 0.05 401, [...], 0.95 0.95 419, 1 1 420),
                   (1 1 1080, 0.95 0.95 1181, [...], 0.05 0.05 1199, 0 0 1200))

instead I get

MULTILINESTRING Z ((0 0 800,0.05 0.05 800),(0.05 0.05 800,0.1 0.1 800),(0.1 0.1 800,0.15 0.15 800),(0.15 0.15 800,0.2 0.2 800),(0.2 0.2 800,0.25 0.25 800),(0.25 0.25 800,0.3 0.3 800),(0.3 0.3 800,0.35 0.35 800),(0.35 0.35 800,0.4 0.4 800),(0.4 0.4 800,0.45 (...)

[pramsey]

I can see how it's possible that multilinestring confuses things… could you try and create a simpler (smaller) example using multilinestrings?

[TBL]

You mean something like this?

Original query with a single line:

with linestring as (
  select st_makeline(
    st_geomfromtext('LINESTRING Z (-1 -1 0,-0.5 -0.5 1,0 0 2,0.5 0.5 3,1 1 4,1.5 1.5 5)', 4326),
    st_geomfromtext('LINESTRING Z (2 2 6,1.5 1.5 7,1 1 8,0.5 0.5 9,0 0 10)', 4326)
  ) as geom
),
polygon as (select ST_GeomFromText('POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))', 4326) as geom)
select st_astext(st_intersection(linestring.geom, polygon.geom)) from linestring, polygon

Result

"MULTILINESTRING Z ((0 0 6,0.5 0.5 6),(0.5 0.5 6,1 1 6))"

Query with a multilinestring:

with linestring as (
  select st_collect(
    st_geomfromtext('LINESTRING Z (-1 -1 0,-0.5 -0.5 1,0 0 2,0.5 0.5 3,1 1 4,1.5 1.5 5)', 4326),
    st_geomfromtext('LINESTRING Z (2 2 6,1.5 1.5 7,1 1 8,0.5 0.5 9,0 0 10)', 4326)
  ) as geom
),
polygon as (select ST_GeomFromText('POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))', 4326) as geom)
select st_astext(st_intersection(linestring.geom, polygon.geom)) from linestring, polygon

Result as above:

"MULTILINESTRING Z ((0 0 6,0.5 0.5 6),(0.5 0.5 6,1 1 6))"

And with two distinct linestrings:

with linestring as (
  select
    (select st_geomfromtext('LINESTRING Z (-1 -1 0,-0.5 -0.5 1,0 0 2,0.5 0.5 3,1 1 4,1.5 1.5 5)', 4326)) as geom
  union all
    (select st_geomfromtext('LINESTRING Z (2 2 6,1.5 1.5 7,1 1 8,0.5 0.5 9,0 0 10)', 4326))
),
polygon as (select ST_GeomFromText('POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))', 4326) as geom)
select st_astext(st_intersection(linestring.geom, polygon.geom)) from linestring, polygon

Result is as expected:

"LINESTRING Z (0 0 2,0.5 0.5 3,1 1 4)"
"LINESTRING Z (1 1 8,0.5 0.5 9,0 0 10)"

[pramsey]

Looks like everything is working as expected… the intersection code has to figure out what to do in the cases where an output coordinate has two *different* z values (which is what happens in both your input examples). The answer, longstanding, is that it takes the average of the input Z values, which is what we're seeing in your examples.

[TBL]

I could imagine something like that when I input a unordered set of points. But my input is a linestring which is a ordered list of points. IMHO ST_Intersection shouldn't try to aggregate anything.

And even if this is the desired behaviour it should be documented.

[TBL]

And BTW, the following example has no common 2d points and the output is still wrong…:

with linestring as (
  select st_collect(
    st_geomfromtext('LINESTRING Z (-1 -1 0,-0.5 -0.5 1,0 0 2,0.5 0.5 3,1 1 4,1.5 1.5 5)', 4326),
    st_geomfromtext('LINESTRING Z (2.01 2.01 6,1.51 1.51 7,1.01 1.01 8,0.51 0.51 9,0.01 0.01 10)', 4326)
  ) as geom
),
polygon as (select ST_GeomFromText('POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))', 4326) as geom)
select st_astext(st_intersection(linestring.geom, polygon.geom)) from linestring, polygon

[pramsey]

No, sorry, it's still just averaging. Take the value at (0,0) for example.

The first line knows that Z at 0,0 is 2 and at 0.5,0.5 is 3, so interpolate, the value at 0.01,0.01 is 2.02

The second line knows that Z at 0.01,0.01 is 10.

The average of those values is 6.01, which is what we see.

[robe]

I think we have a couple of functions besides ST_Intersection where similar behavior exists from recollection:

ST_ConvexHull, ST_Union

am I missing any.

[robe]

I think this is sorta already resulting in expected behavior if you have sfcgal support.

Just for example:

set postgis.backend = sfcgal;
with linestring as (
  select st_makeline(
    (select st_setsrid(st_makeline(st_makepoint((a-400)/20.0,(a-400)/20.0,a)),4326) from generate_series(0,800) as a),
    (select st_setsrid(st_makeline(st_makepoint(-(a-1200)/20.0,-(a-1200)/20.0,a)),4326) from generate_series(801,1600) as a)
  ) as geom
),
polygon as (select ST_GeomFromText('POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))', 4326) as geom)
select st_astext(st_intersection(linestring.geom, polygon.geom)) from linestring, polygon

output is 2D as many people would expect from a 2D function:

MULTILINESTRING((0.05 0.05,0 0),(0.1 0.1,0.05 0.05),(0.15 0.15,0.1 0.1),(0.2 0.2,0.15 0.15),(0.25 0.25,0.2 0.2),(0.3 0.3,0.25 0.25),(0.35 0.35,0.3 0.3),(0.4 0.4,0.35 0.35),(0.45 0.45,0.4 0.4),(0.5 0.5,0.45 0.45),(0.55 0.55,0.5 0.5),(0.6 0.6,0.55 0.55),(0.65 0.65,0.6 0.6),(0.7 0.7,0.65 0.65),(0.75 0.75,0.7 0.7),(0.8 0.8,0.75 0.75),(0.85 0.85,0.8 0.8),(0.9 0.9,0.85 0.85),(0.95 0.95,0.9 0.9),(1 1,0.95 0.95))

now if I use the ST_3DIntersection

set postgis.backend=sfcgal
with linestring as (
  select st_makeline(
    (select st_setsrid(st_makeline(st_makepoint((a-400)/20.0,(a-400)/20.0,a)),4326) from generate_series(0,800) as a),
    (select st_setsrid(st_makeline(st_makepoint(-(a-1200)/20.0,-(a-1200)/20.0,a)),4326) from generate_series(801,1600) as a)
  ) as geom
),
polygon as (select ST_GeomFromText('POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))', 4326) as geom)
select st_astext(st_3dintersection(linestring.geom, polygon.geom)), ST_3DIntersects(linestring.geom, polygon.geom) from linestring, polygon

I get

       st_astext         | st_3dintersects
-------------------------+-----------------
GEOMETRYCOLLECTION EMPTY | f

I tried swapping out change

set postgis.backend = geos;

and both the built in native ST_3DIntersects agrees that these don't intersect.

Now if I take a slightly modified that should intersect it does

with linestring as (
  select st_makeline(
    st_geomfromtext('LINESTRING Z (-1 -1 0,-0.5 -0.5 1,0 0 2,0.5 0.5 3,1 1 4,1.5 1.5 5)', 4326),
    st_geomfromtext('LINESTRING Z (2 2 6,1.5 1.5 7,1 1 8,0.5 0.5 8,0 0 10)', 4326)
  ) as geom
),
polygon as (select ST_GeomFromText('POLYGON((0 0 8, 0 1 8, 1 1 8, 1 0 8, 0 0 8))', 4326) as geom)
select st_astext(st_3dintersection(linestring.geom, polygon.geom)) from linestring, polygon;

— output is —

LINESTRING Z (1 1 8,0.5 0.5 8)

[robe]

I have put a caution on the ST_Intersection function that is does some sort of averaging of Z coordinate and also added sfcgal behavior (for both ST_Intersection and ST_3DIntersection) for comparison. Done at r12753.

I think it needs more work so I'll leave this ticket open.

Once debbie is done building doc, you can see changes here:

​http://postgis.net/docs/manual-dev/ST_Intersection.html

and here: ​http://postgis.net/docs/manual-dev/ST_3DIntersection.html

[robe]

No replied to this closing. Fix was put in trunk (not backported)