2

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.

5

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)
3

The easiest way I can think of to get a correct answer is to use the point that you have found with ST_Intersection, build a vertical line through that and find where that line intersects with the line you want the z-value from.

This vertical line have to be long enough to intersect with your line of interest. You can build it with something like:

SELECT 
ST_MakeLine(ST_MakePoint(ST_X(p), ST_Y(p),0),ST_MakePoint(ST_X(p), ST_Y(p), 100000))
  FROM 
(
    SELECT ST_Intersection(a,b) p, b FROM 
    (
        SELECT 'LINESTRING(2 0 0, 2 6 50)'::geometry a,'LINESTRING(0 3 0, 5 3 0)'::geometry b
    ) lines
) foo

Then to find where this vertical line intersects with geometry a in my example you can use ST_3DClosestPoint. And to get the Z-value you can use ST_Z on the resulting point.

All together will look something like this frommy example:

SELECT 
ST_Z(ST_3dClosestPoint(
a, ST_MakeLine(ST_MakePoint(ST_X(p), ST_Y(p),0),ST_MakePoint(ST_X(p), ST_Y(p), 100000))))
  FROM 
(
    SELECT ST_Intersection(a,b) AS p, a FROM 
    (
        SELECT 'LINESTRING(2 0 0, 2 6 50)'::geometry a,'LINESTRING(0 3 0, 5 3 0)'::geometry b
    ) lines
) foo

Note that just finding the point where the two lines are closest two eachother won't give a correct answer. You can try it by looking at the answer from

SELECT ST_ASEWKT(ST_3dShortestLine(a,b))
  from 
(SELECT 'LINESTRING(2 0 0, 2 6 50)'::geometry a,'LINESTRING(0 3 0, 5 3 0)'::geometry b) lines

Since the lines don't intersect in 3D the shortest line between the two linestrings will not be vertical in most cases and you will get the wrong answer.

0

You could use ST_Intersection(linework.geom, align.geom) as point with ST_LineLocatePoint(interesting_geom, point) as fraction_of_ig and use it get z/m value using ST_LineInterpolatePoint(interesting_geom, fraction_of_ig ).

  • This will not give the right answer in all cases. example is the linestrings: LINESTRING(2 0 0, 2 6 50) and LINESTRING(0 3 0, 5 3 0) – Nicklas Avén Nov 11 '13 at 20:47
  • ST_AsText(ST_Intersection(ST_GeomFromText('LineString(2 0 0, 2 6 50)'),ST_GeomFromText('LineString(0 3 0, 5 3 0)'))) returns "POINT Z (2 3 12.5)".. i didn know that. – simplexio Nov 12 '13 at 15:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.