[GIS] ST_Intersection point of two 3D linestrings – how is the Z value calculated

intersectionpostgis

I have two 3D linestrings (geometry type LINESTRING Z) in a PostGIS database that intersect in 2D space, but not in 3D space. Given that there is technically no 'shared portion' of geometry in 3D, how is the Z value returned by ST_Intersection calculated?

Here is the query I am using:

SELECT feature, str1, (ST_Dump(geom)).geom.ST_Z as z_value
FROM
(SELECT ST_Intersection(linework.geom, align.geom) as geom, linework.feature, linework.str1
FROM align
INNER JOIN linework ON ST_Intersects(align.geom, linework.geom)
WHERE ST_isvalid(align.geom)='t' AND ST_isvalid(linework.geom)='t') q1;

What I would like is a query that returns the z value of the point of intersection but on one of the linestrings only.

linework contains surveyed features that are crossed by an alignment in align – I would like to 'lift' elevations off the surveyed features at the points of intersection but I'm not sure how to go about it.

Best Answer

For ST_Intersection, PostGIS uses one of two libraries. Here is the source code for that.

  • GEOS is by far the most common library used by PostGIS for this operation. GEOS is similar to JTS, in that it is a 2D library, with limited interpretations of the Z dimension sort-of tacked on top.
  • SFCGAL is Simple Features wrapper to CGAL, and is a recent development introduced with PostGIS 2.1. This geometry library is more 3D-aware, so results may differ than from GEOS.

As I only have the GEOS version, I can only say what this library is doing. Take two linestrings that intersect in 2D:

SELECT ST_AsText(ST_Intersection('LINESTRING Z (0 0 0, 2 2 2)',
                                 'LINESTRING Z (0 1 4, 1 0 4)'));
       st_astext
------------------------
 POINT Z (0.5 0.5 2.25)
(1 row)

The intersection point for each has Z values of 0.5 (first quarter length) and 4.0 (it's a flat line at 4.0). The average of these two numbers is 2.25. So we can say it first determines the 2D location of the intersection, does a linear interpolation of the Z value for each linestring, then takes the average of the two Z values.

To interpolate the Z value from only one of the linestrings, use the linear referencing functions ST_LineInterpolatePoint and ST_LineLocatePoint (for PostGIS 2.0 and before, these are called ST_Line_Interpolate_Point and ST_Line_Locate_Point).

SELECT ST_AsText(ST_LineInterpolatePoint(A, ST_LineLocatePoint(A, ST_Intersection(A, B))))
FROM (SELECT 'LINESTRING Z (0 0 0, 2 2 2)'::geometry AS A,
             'LINESTRING Z (0 1 4, 1 0 4)'::geometry AS B) AS f;
       st_astext
-----------------------
 POINT Z (0.5 0.5 0.5)

To simplify things, here's a function that will return an interpolated point using the Z value from the first linestring:

CREATE FUNCTION ST_IntersectsFirstZ(linestringA geometry, linestringB geometry)
RETURNS geometry AS
'SELECT ST_Line_Interpolate_Point($1, ST_Line_Locate_Point($1, ST_Intersection($1, $2)))'
LANGUAGE sql IMMUTABLE;

And using it:

SELECT ST_AsText(ST_IntersectsFirstZ(A, B)) AS AB,
       ST_AsText(ST_IntersectsFirstZ(B, A)) AS BA
FROM (SELECT 'LINESTRING Z (0 0 0, 2 2 2)'::geometry AS A,
             'LINESTRING Z (0 1 4, 1 0 4)'::geometry AS B) AS f;
-[ RECORD 1 ]-------------
ab | POINT Z (0.5 0.5 0.5)
ba | POINT Z (0.5 0.5 4)
Related Question