Query Z at intersecting point locations from a MultiPolygonZ?

Question

I create a MultiPolygonZ geometry with ST_DelaunayTriangles().

select st_setsrid(
           ST_CollectionExtract(
               ST_DelaunayTriangles(
                   ST_GeomFromText('MULTIPOINT Z(-1.293574 53.478486 1293, -1.077978 53.734323 2527, -1.077978 53.221097 2250, -0.862383 53.478486 2250, -1.50917 53.734323 2595, -1.50917 53.221097 3277, -1.293574 53.988613 3351, -1.724765 53.478486 3813, -1.293574 52.962151 3868, -0.862383 53.988613 3978, -0.646787 53.734323 3202, -0.646787 53.221097 3134, -0.862383 52.962151 3794, -1.724765 53.988613 4598, -1.077978 54.241359 5150, -1.724765 52.962151 4834, -1.940361 53.734323 4437, -1.940361 53.221097 5463, -0.431191 53.478486 3761, -1.50917 54.241359 4178, -1.077978 52.701646 4822, -1.50917 52.701646 4646, -0.431191 53.988613 4824, -0.646787 54.241359 5346, -2.155957 53.478486 5592, -1.293574 53.563937 1337, -0.862383 53.563937 2057, -1.077978 53.307066 1806, -0.970181 53.435696 2009, -1.185776 53.435696 1244, -1.401372 53.435696 2073, -0.754585 53.435696 3145, -0.970181 53.178047 2517, -1.185776 53.178047 2894, -1.293574 53.307066 1863, -0.862383 53.307066 2434, -1.077978 53.04864 3172, -1.077978 53.563937 1455, -1.50917 53.563937 2635, -0.646787 53.563937 2942, -0.970181 53.691791 2027, -1.185776 53.691791 1905, -1.401372 53.691791 2084, -0.754585 53.691791 3251, -1.50917 53.307066 2839, -0.646787 53.307066 3377, -0.862383 53.04864 3200, -1.293574 53.04864 2909, -1.401372 53.178047 2711, -0.754585 53.178047 3047, -1.401372 53.734323 2364, -1.185776 53.734323 2166, -0.970181 53.734323 2515, -1.293574 53.606598 1302, -0.862383 53.606598 2263, -1.077978 53.606598 1831, -1.401372 53.478486 1634, -1.185776 53.478486 854, -0.970181 53.478486 1832, -1.293574 53.349986 1634, -1.077978 53.349986 1542, -1.024079 53.670509 2299, -1.131877 53.670509 1863, -1.455271 53.670509 2284, -1.347473 53.670509 1961, -1.239675 53.670509 1539, -0.916282 53.670509 2431, -1.024079 53.542591 1247, -0.916282 53.542591 2307, -1.239675 53.542591 908, -1.347473 53.542591 1365, -0.808484 53.542591 2264, -1.131877 53.542591 1210, -1.024079 53.414284 1496, -1.131877 53.414284 1627, -1.239675 53.691791 2045, -1.293574 53.627913 1519, -1.347473 53.563937 1537, -1.239675 53.563937 1037, -1.131877 53.563937 1433, -1.024079 53.563937 1611, -1.401372 53.499865 1824, -1.185776 53.499865 539, -1.131877 53.435696 1555, -1.024079 53.435696 1632, -1.293574 53.371429 1846, -1.077978 53.371429 1583, -1.212726 53.659864 1708, -1.266625 53.659864 1325, -1.266625 53.595937 1265, -1.320523 53.595937 1624, -0.99713 53.531913 1928, -1.104928 53.531913 1280, -1.212726 53.531913 895, -1.374422 53.531913 1575, -1.320523 53.531913 1195, -1.266625 53.531913 860, -1.158827 53.531913 919, -1.051029 53.531913 1651, -1.158827 53.467792 1175)')
               ),
           3),
        4326)::geometry(multipolygonz, 4326) geom

I use st_dwithin() to find point geometries in a Postgis table.

Is it possible to assign the Z value at the intersections to the points without going to raster first?

I tried st_intersection() like so:

SELECT st_intersection(eu_hx_8k.geomcntr, delaunay.geom) geom
FROM
eu_hx_8k,
(subquery) delaunay
WHERE st_dwithin(delaunay.geom, eu_hx_8k.geomcntr, 0);

This finds the intersecting points but st_intersection() fails with a TopologyException.

[XX000] ERROR: Error performing intersection: TopologyException: Input geom 1 is invalid: Self-intersection at or near point -2.1559569999999999 53.478485999999997 5592 at -2.1559569999999999 53.478485999999997 5592

Edit:

I am getting somewhere by breaking the MultiPolygonZ into PolygonZ using ST_Dump.

I then intersect PointZ which I create with ST_MakePoint().

However, only the PointZ which intersect the boundary of the PolygonZ return a Z.

SELECT
  eu_hx_8k.geomcntr,
  ST_Z(ST_Intersection(ST_SetSRID(ST_MakePoint(eu_hx_8k.lon, eu_hx_8k.lat, 0), 4326), delaunay.geom)) z
FROM
  eu_hx_8k,
  (SELECT (ST_Dump(
      ST_SetSRID(
          ST_CollectionExtract(
              ST_DelaunayTriangles(
                  ST_GeomFromText(
                      'MULTIPOINT Z(-1.293574 53.478486 1293, -1.077978 53.734323 2527, -1.077978 53.221097 2250, -0.862383 53.478486 2250, -1.50917 53.734323 2595, -1.50917 53.221097 3277, -1.293574 53.988613 3351, -1.724765 53.478486 3813, -1.293574 52.962151 3868, -0.862383 53.988613 3978, -0.646787 53.734323 3202, -0.646787 53.221097 3134, -0.862383 52.962151 3794, -1.724765 53.988613 4598, -1.077978 54.241359 5150, -1.724765 52.962151 4834, -1.940361 53.734323 4437, -1.940361 53.221097 5463, -0.431191 53.478486 3761, -1.50917 54.241359 4178, -1.077978 52.701646 4822, -1.50917 52.701646 4646, -0.431191 53.988613 4824, -0.646787 54.241359 5346, -2.155957 53.478486 5592, -1.293574 53.563937 1337, -0.862383 53.563937 2057, -1.077978 53.307066 1806, -0.970181 53.435696 2009, -1.185776 53.435696 1244, -1.401372 53.435696 2073, -0.754585 53.435696 3145, -0.970181 53.178047 2517, -1.185776 53.178047 2894, -1.293574 53.307066 1863, -0.862383 53.307066 2434, -1.077978 53.04864 3172, -1.077978 53.563937 1455, -1.50917 53.563937 2635, -0.646787 53.563937 2942, -0.970181 53.691791 2027, -1.185776 53.691791 1905, -1.401372 53.691791 2084, -0.754585 53.691791 3251, -1.50917 53.307066 2839, -0.646787 53.307066 3377, -0.862383 53.04864 3200, -1.293574 53.04864 2909, -1.401372 53.178047 2711, -0.754585 53.178047 3047, -1.401372 53.734323 2364, -1.185776 53.734323 2166, -0.970181 53.734323 2515, -1.293574 53.606598 1302, -0.862383 53.606598 2263, -1.077978 53.606598 1831, -1.401372 53.478486 1634, -1.185776 53.478486 854, -0.970181 53.478486 1832, -1.293574 53.349986 1634, -1.077978 53.349986 1542, -1.024079 53.670509 2299, -1.131877 53.670509 1863, -1.455271 53.670509 2284, -1.347473 53.670509 1961, -1.239675 53.670509 1539, -0.916282 53.670509 2431, -1.024079 53.542591 1247, -0.916282 53.542591 2307, -1.239675 53.542591 908, -1.347473 53.542591 1365, -0.808484 53.542591 2264, -1.131877 53.542591 1210, -1.024079 53.414284 1496, -1.131877 53.414284 1627, -1.239675 53.691791 2045, -1.293574 53.627913 1519, -1.347473 53.563937 1537, -1.239675 53.563937 1037, -1.131877 53.563937 1433, -1.024079 53.563937 1611, -1.401372 53.499865 1824, -1.185776 53.499865 539, -1.131877 53.435696 1555, -1.024079 53.435696 1632, -1.293574 53.371429 1846, -1.077978 53.371429 1583, -1.212726 53.659864 1708, -1.266625 53.659864 1325, -1.266625 53.595937 1265, -1.320523 53.595937 1624, -0.99713 53.531913 1928, -1.104928 53.531913 1280, -1.212726 53.531913 895, -1.374422 53.531913 1575, -1.320523 53.531913 1195, -1.266625 53.531913 860, -1.158827 53.531913 919, -1.051029 53.531913 1651, -1.158827 53.467792 1175)')
              ),
              3),
          4326))).geom :: GEOMETRY(polygonz, 4326) geom) delaunay
WHERE st_dwithin(delaunay.geom, eu_hx_8k.geomcntr, 0);

Edit:

I was able to solve this with turf in node.js. I still want to keep this question open as there might be a PostGIS solution which is more performant but then again TurfJS itself is lightning quick solving this problem.

Answer 1

A colleague send me this bit of Postgres SQLin order to built a function which would do the same calculation directly in the database. Unfortunately I won't have time now to test whether this approach can be faster than TurfJS. I just want to put this here for reference. Perhaps somebody else will find this helpful and writes a PostGIS function.

declare @pointA geometry; 
declare @pointB geometry; 
declare @pointC geometry; 
declare @pointAweight float; 
declare @pointBweight float; 
declare @pointCweight float; 


set @pointA=geometry::STGeomFromText('Point(175.0305935740471 -39.924665194652604)', 4326)
set @pointB=geometry::STGeomFromText('Point(175.03033608198166 -39.924387504970255)', 4326) 
set @pointC=geometry::STGeomFromText('Point(175.0301563739777 -39.92449035313209)', 4326)
set @pointAweight=5.1
set @pointBweight=3.7
set @pointCweight=2.9

---- interpolate for new point (175.0303 -39.9245)
-- first calculate area of full triangle
declare @ftri geometry;
declare @ftri_area float;
set @ftri=geometry::STGeomFromText('POLYGON((175.0305935740471 -39.924665194652604,175.03033608198166 -39.924387504970255,175.0301563739777 -39.92449035313209,175.0305935740471 -39.924665194652604))', 4326)
set @[email protected]()

-- now find the areas of the triangle cut up into regions (note the omiited point relates to the point to apply the final weighting) 
declare @triA geometry;
declare @triA_area float;
set @triA=geometry::STGeomFromText('POLYGON((175.0303 -39.9245,175.03033608198166 -39.924387504970255,175.0301563739777 -39.92449035313209,175.0303 -39.9245))', 4326)
set @[email protected]()


declare @triB geometry;
declare @triB_area float;
set @triB=geometry::STGeomFromText('POLYGON((175.0305935740471 -39.924665194652604,175.0303 -39.9245,175.0301563739777 -39.92449035313209,175.0305935740471 -39.924665194652604))', 4326)
set @[email protected]()
select @triB

declare @triC geometry;
declare @triC_area float;
set @triC=geometry::STGeomFromText('POLYGON((175.0305935740471 -39.924665194652604,175.03033608198166 -39.924387504970255,175.0303 -39.9245,175.0305935740471 -39.924665194652604))', 4326)
set @[email protected]()

-- build interpolated value

select @pointAweight*@triA_area/@ftri_area+@pointBweight*@triB_area/@ftri_area+@pointCweight*@triC_area/@ftri_area  'interpolated_value'

Answer 2

You can solve this using built-in PostGIS functions if you have enabled the SFCGAL extension.

Using ST_3DIntersection, you can "poke" your polygonZ with vertical lines created from the point locations that you wish to intersect.

Here is a quick example using a simple triangular polygonZ and a 2d point:

WITH a_pgz AS --Input polygonZ
(SELECT ST_GeomFromText('POLYGON Z((2.5 5 10.3,-2 20 25,20 0 15,2.5 5 10.3))') AS geom),
b_point2d AS --Input 2d point
(SELECT ST_GeomFromText('POINT (5.5 6.7)') AS geom),
c_linestringz AS --Build vertical linestringZ from 2d point by adding z values: Floor(ST_ZMin(a_pgz.geom)))-1 = z of linestring start point, Round(ST_ZMax(a_pgz.geom)))+1 = z of linestring end point
(SELECT ST_MakeLine(ST_Translate(ST_Force3D(b_point2d.geom), 0, 0, Floor(ST_ZMin(a_pgz.geom))-1), ST_Translate(ST_Force3D(b_point2d.geom), 0, 0, Round(ST_ZMax(a_pgz.geom))+1)) AS geom
FROM a_pgz, b_point2d)
--Get the Z value of the intersection of the polygonZ and the linestringZ
SELECT ST_Z(ST_3dIntersection(a_pgz.geom, c_linestringz.geom)) AS z_pt_on_pgz
FROM a_pgz, c_linestringz