postgresql - PostGIS -> LineString + PointZ = LineStringZ

DATA

I have 2 tables in PostGIS:

• linestable with MultiLineString geom


• surveypointtable with Point + elevation attribute

GOAL

The goal is to rebuild the lines as 3D replacing line's vertices with the closets surveypoint.

LOGIC

1. 
   

   Rebuild points as PointZ and get the collection of all points

   
   

   ST_Collect(ST_MakePoint(ST_X(geom),ST_Y(geom),elevation_att))
   



2. 
   

   Get line's vertices as points using

   
   

   ST_DumpPoints(geom).geom
   




3. 
   

   Find closest surveypoint using

   
   

   ST_ClosestPoint(surveypoints.geom, linepoints.geom)
   



4. 
   

   Check if the closestpoint is within 1m using

   
   

   ST_DWithin(surveypoints.geom, linepoints.geom, 1)
   



5. 
   

   Build line using only the closets surveypoints that are within 1m tolerance.

   
   
   

   ST_MakeLine(ST_Collect(closestpoints.geom))
   

PROBLEM

When I use ST_ClosestPoint(p.3) it drops the z-index of the point and as a result I end up with a 2D line.

When I use ST_3DClosestPoint it takes wrong points - i.e. if I had point A 1m away in 2D from linepoint but at the same elevation and point B at the exact same coordinates as linepoint but with 10m difference in elevation it would take point A.

CODE

    DROP TABLE IF EXISTS newlines;



    CREATE TABLE newlines AS
    SELECT gid,
          ST_SetSRID(ST_MakeLine(geom),27700)
    FROM
      (SELECT gid,
              path,
              ST_AsText(
                ST_Collect(
                  closestsurveypoint)) AS geom
      FROM

        ( SELECT linepoints.gid AS gid,
                  linepoints.path AS path,
                  linepoints.geom AS linepoint,
                  (CASE ST_DWithin(
                        ST_SetSRID(surveypoints.geom, 27700), 
                        ST_SetSRID(linepoints.geom, 27700), 
                        1)
                      WHEN TRUE THEN ST_AsText(
                        ST_ClosestPoint(
                          ST_SetSRID(surveypoints.geom, 27700), 

                          ST_SetSRID(linepoints.geom, 27700) ))
                  END) AS closestsurveypoint
          FROM
            ( SELECT gid,
                    color, 
                    (ST_DumpPoints(geom)).path AS path,
                    ST_AsText((ST_DumpPoints(geom)).geom) AS geom
            FROM linestable) AS linepoints,

            ( SELECT ST_AsText(

              ST_Collect(
                ST_MakePoint(
                  ST_X(geom), 
                  ST_Y(geom), 
                  elevation))) AS geom
            FROM surveypointstable) AS surveypoints ) AS newlinepoints
      GROUP BY gid,
                path
      ORDER BY gid,
                path) AS newlines

    GROUP BY gid

Answer

Don't create POINTZ to extrapolate the elevation; instead, find the nearest neighbor of your surveypointtable to each vertex of linestable first, then create POINTZ from each vertex with that respective elevation attribute and stitch the lines back together via row number and the path sequence of ST_DumpPoints:

WITH
    seqs AS (
        SELECT --pts., pts.,
               mln_seq,
               pts.path[1] AS elem_seq,
               pts.path[2] AS pnt_seq,
               pts.geom AS geom

        FROM (
            SELECT --, , 
                   ROW_NUMBER() OVER() AS mln_seq,
                   (ST_DumpPoints(geom)).*
            FROM linestable
        ) AS pts
    )

SELECT DISTINCT
       --s., s.,

       ST_MakeLine(ST_MakePoint(ST_X(s.geom), ST_Y(s.geom), q.elevation)) OVER(PARTITION BY s.mln_seq, s.elem_seq ORDER BY s.elem_seq) AS geom
FROM seqs AS s
JOIN LATERAL (
    SELECT elevation
    FROM surveypointtable
    ORDER BY geom <-> s.geom
    LIMIT 1
) AS q
ON true;

I included , for the linestable; replace and uncomment if you intent to keep attributes from your original lines.

Note that you will have to set SRIDs as necessary; have a look in the tags nearest-neighbor and knn for postgis to find info on the used structure (LATERAL JOIN), including detailed description of proper CRS/unit usage and the <-> KNN operator.

Note also that I assume your initial lines are MULTILINESTRINGs, since you said so; the result will be simple LINESTRINGZs, as your title suggests. Honestly, I am not sure what you have and expect...but this should work with LINESTIRNGs as input if you remove all references to pnt_seq from above (i.e. pts.path[2], as there will only be one sequence number in path), and you can ST_Collect the result grouped by mln_seq if you need MULTILINESTRINGZs; then, however, you cannot keep attributes as simple as above.