2

In a nutshell

I'd like to compute the intersection of a 3D plane (defined by 3 points) and a terrain (given as a raster DEM GeoTIFF), and more particularly the part between two lines of that plane.

Inputs

I have a (GeoTIFF) raster DEM representing terrain height in a given projected (Cartesian) coordinate system.

Here's an hillshade representation of this DEM: Hillshade of the DEM with three point

I have 3 points in space above this DEM representing an arbitrarily oriented plane. Each point has a known (x,y,z) position in the same coordinate system as the terrain. (These points are currently stored in a PostGIS database as POINT geometries).

I want to compute the intersection of this plane with the terrain, especially the part (depicted in cyan in the image) between two of the three edges (depicted by yellow lines) of the triangle defined by these 3 points:

Plane intersection with the terrain

Upon some thoughts, it may not be so trivial.

What I'm 90% sure (based on intuitions) is that:

  • I probably need to compute triangles from the raster DEM
  • It's only the intersections (dark blue dots on the image below) between the plane and all the edges of these terrain triangles that is of interest (lines between these point can be interpolated; in cyan), as explained by this closeup on terrain triangles:

terrain triangle intersections zoomed

  • It's probably totally useless to test each and every edge of the terrain triangles, over the whole DEM, as only a tiny set of these will actually intersect, but if you don't know them in advance, you probably have to establish kind of a (tricky?) optimized search procedure to quickly excludes the majority of them
  • The computation time will certainly be proportional to the size of the DEM
  • All these points make me think that I probably need a low-level implementation (languages of the C family or similar, a pure PostGIS implementation would be truly amazing) but an implementation that is "user friendly" (to be used by a non C-expert)

I'm 100% sure that:

  • I need to run this computation on a headless Ubuntu machine

Directions I explored:

I've searched QGIS plugins deeply, but didn't found anything relevant enough to solve this.

I've also searched PostGIS functions without much success, except perhaps on the side of the sfcgal extension but I was not able to figure out a way to do such computation, especially when the terrain is stored as a file. And on the other hand, as the plane is not "outlined", standard functions such as ST_Intersection could not help much.

I also looked a some viewshed algorithms (GDAL, SAGA), hoping to find some that could use a 3D line to delineate the result, without much success either.

And I lost myself in Python libraries...

I want to first turn the 90% of assumptions above into 100% certainty, and then find a friendly way to go? But which one and how; have you ever successfully managed such a calculation without implementing all the steps "by hand" (which in my case will certainly takes ~hours to run if I try in pure Python)?

Improve this question
4
  • 1
    I'd definitely try to jump out of the world of GIS applications and into common 3D tools as fast as possible (to be as fast and efficient as possible). If you want to stay with Python, check trimesh or pymesh. For triangulation try tin-terrain or pydelatin.
    – bugmenot123
    Commented Jul 6, 2021 at 19:36
  • 2
    Very basic task if you use rasters instead. Compute new surface and find where dem is above it.
    – FelixIP
    Commented Jul 6, 2021 at 20:20
  • @FelixIP - I was thinking along these lines as well, but is there a simple way to create a raster corresponding to the sloping plane?
    – Llaves
    Commented Jul 6, 2021 at 22:14
  • Thanks bugmenot123, I'll explore these promising possibilities. And thanks also to FelixP for this seducing idea, but which has one major drawback; I do have lots of vertical to "close-to-vertical" planes. This will inevitably lead to instabilities in the numerical calculations.
    – swiss_knight
    Commented Jul 7, 2021 at 17:05

1 Answer 1

Reset to default
1

Picture below shows:

  • 3 points labelled with 'elevation'
  • contours of original DEM
  • contours of plain surface made from 3 lines enter image description here

Points coordinates:

enter image description here

Found equation of a surface using regression analysis in Excel:

enter image description here

The rest was very ArcGIS specific:

  • compute raster of X
  • compute raster of Y

Aplly above equation in raster calculator to produce plain surface raster and use calculator to find areas where it sits below original.

You can apply it to all sort of things, e.g. design of forestry roads...

Improve this answer
3
  • 1
    Could you please elaborate with more detail? Right now with your answer it is not possible to follow your steps and concepts to understand if it solves the issue. What exactly you do in Excel? What exactly do you do in ArcGIS?
    – bugmenot123
    Commented Jul 7, 2021 at 15:20
  • As you are using a regression, I think your coefficients are not as correct as there would be when solving the actual plane equation (I do have small differences on my side). You definitely don't need a regression here. And using Excel is doubly not compatible with a headless Linux server. But I got the idea, thanks a lot!
    – swiss_knight
    Commented Jul 7, 2021 at 17:09
  • It IS about idea, mate. Implementation is up to you. I use Excel to solve a set of 3 linear equation with 3 unknowns because I won't bother to spend more time on answer.
    – FelixIP
    Commented Jul 7, 2021 at 19:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.