Tuesday, 4 April 2017

postgis - ST_Intersection point of two 3D linestrings - how is the Z value calculated?


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

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...