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.
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)
No comments:
Post a Comment