[GIS] 3D Points within 3D Polygon in PostGIS

3dpointpolygonpostgisvolume

I've been working with PostGIS for a while now, and I'm able to run queries over massive datasets reasonably easy and fast. However, I do have an issue now that I don't know how to solve. I did find some workarounds, although they mean lower performance, higher execution times, and more complex queries.

The objective is straightforward: to get all the points within a volume.
But, the volume is not a simple extrusion of a polygon, it is instead of the difference of a higher volume minus some smaller volumes here and there.

My first attempt: to tackle the problem was inspired by this question, although it didn't work:

  1. Create the master polygon, using a third axis (Z) as the base for the volume:

ST_Polygon(ST_GeomFromText('LINESTRING(Longs, Lats, Zs...)'), 4326);

  1. Extrude the polygon using the height of the volume, getting a POLYHEDRALSURFACE as a resulting geometry:

ST_Extrude(Previous_polygon, x=0, y=0, z=height);

  1. Having a table with the information regarding each master volume and its small volumes to be removed from. I created a new geometry removing all these small volumes from the master volume:

ST_3DDifference(Master volume, combination of smaller volumes to be removed);

  1. Query the table with the positions and joining the resulting master volume on

3dpoints.geometry &&& Master_Volume.geometry.

As I said, it didn't work because the points were also found outside the master volume. Not sure why! I guess I'll try to come back to that option should the following solution be too much of a hassle.

The second attempt I tried was:

  1. Query the table with the positions straightaway, merging the volume as a 2D polygon, and including the filter on the Z-axis, something like:

ST_Contains(Master_Polygon.geometry, 2dpoints.geometry) AND 2dpoints.Z BETWEEN Master_Polygon.MinZ AND Master_Polygon.MaxZ

  1. Of course, this included all the point within the smaller volumes I didn't want to consider, so then I added the condition to remove those points:

ST_Contains(Master_Polygon.geometry, 2dpoints.geometry) AND 2dpoints.Z BETWEEN Master_Polygon.MinZ AND Master_Polygon.MaxZ AND NOT (ST_Contains(Smaller_Polygon.geometry, 2dpoints.geometry) AND 2dpoints.Z BETWEEN Smaller_Polygon.MinZ AND Smaller_Polygon.MaxZ)

And it does work, although only when the smaller volume (polygon) is only one. Should the lower volumes be more than one, then the condition is not enough as it would create a different set of points based on each condition. This gives all the points within the master volume but not in the lower volume A would be returned or all the points within the master volume but not in the lower volume B.

Question:
Is there a way to get all the points in the master volume by removing all those points with a condition which are in any of the smaller volumes?

I apologize for the long post, but its due to the complicated explanation of the issue.
I did find a solution, and it's by splitting the master volume into smaller volumes and merging it after the query. But it is not quite right, especially when I will have to update the geometries and make the master volume to change accordingly automatically.

Best Answer

I played around with some test setups a while ago; it seems ST_3DUnion, ST_3DDifference and ST_3DIntersection work on Volumes; you could try sth. like this:

WITH
  volume AS (
    SELECT ST_3DDifference(
              <master_volume>, 
              ST_3DUnion(<smaller_volumes>)
           ) AS geom
  )
SELECT ST_3DIntersection(vol.geom, pts.geom) AS geom
FROM volume AS vol,
     <points> AS pts

Since the spatial relation filters won´t work, a selective JOIN is impossible; thus this will take damn long and returns GEOMETRYCOLLECTION EMPTY for all non-intersecting points that you then need to filter out.

Just wanted to share; I´m not sure if this is actually going to work, maybe someone would know exactly if and why...

It seems to return correct results with simple test setups like

SELECT ST_AsText(
            ST_3DIntersection(
                ST_3DDifference(
                    ST_MakeSolid(ST_Extrude(ST_Buffer('POINTZ(0 0 0)'::geometry, 20), 0, 0, 20)),
                    ST_3DUnion(ST_MakeSolid(ST_Extrude(ST_Buffer('POINTZ(0 0 0)'::geometry, 1), 0, 0, 10)),
                               ST_MakeSolid(ST_Extrude(ST_Buffer(ST_GeomFromText('POINTZ(0 0 5)'), 1), 0, 0, 10)))
                ),
                'POINTZ(0 0 5)'::geometry     --'POINTZ(2 0 5)'::geometry
            )
       ) AS geom

but that could also just be completely random...

Related Question