3

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
  • Have you tried using ST_Force3D? I have found with certain 3D operations, I had to wrap everything in such calls to get correct behaviour. – John Powell Oct 18 '18 at 10:30
4

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.<some_column>, pts.<some_other_column>,
               mln_seq,
               pts.path[1] AS elem_seq,
               pts.path[2] AS pnt_seq,
               pts.geom AS geom
        FROM (
            SELECT --<some_column>, <some_other_column>, 
                   ROW_NUMBER() OVER() AS mln_seq,
                   (ST_DumpPoints(geom)).*
            FROM linestable
        ) AS pts
    )

SELECT DISTINCT
       --s.<some_column>, s.<some_other_column>,
       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 <some_column>, <some_other_column> 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 and for 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.

  • ah, btw. if you want to create a new table, you can simply add CREATE TABLE AS ... before the SELECT DISTINCT ..., i.e. you can use CTEs for table creation queries. – ThingumaBob Oct 18 '18 at 11:39
  • WoW! This is BRILLIANT! It works like a charm, thank you very much. – Łukasz Kłosiński Oct 18 '18 at 11:48
  • I have one more question - what happens if it doesn't find a point or point is with NULL elevation? – Łukasz Kłosiński Oct 18 '18 at 11:49
  • And answering your question - I have MultiLineString in my linestable but can work with LineStrings as the query's output. – Łukasz Kłosiński Oct 18 '18 at 11:54
  • 1
    @ŁukaszKłosiński it will always find a point...if there is at least one point in your surveypointtable; you even have the option to define a WHERE filter in the LATERAL sub-query, e.g. if you want to consider only points with a matching column value etc. for each line. as to the NULL value: it will probably throw an error in ST_MakePoint; use SELECT COALESCE(elevation, <default_value>) ... instead of SELECT elevation ... in the sub-query to catch that case. – ThingumaBob Oct 18 '18 at 11:59

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.