3

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.

enter image description here

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.

enter image description here

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.

enter image description here

  • It is difficult to understand what you are trying to do. Can you elaborate further on what is your use case at hand and what would be your desired result? – tilt Apr 20 '17 at 19:46
  • The desired result is what is called a planepoint in turf.js turfjs.org/docs/#planepoint . I will actually try use this since I am having a node.js backend and turf.js is usually super fast. – Dennis Bauszus Apr 26 '17 at 13:30
0

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 @ftri_area=@ftri.STArea()

-- 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 @triA_area=@triA.STArea()


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 @triB_area=@triB.STArea()
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 @triC_area=@triC.STArea()

-- build interpolated value

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

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

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.