extend rasterize_terrain()

[AdrianKriger]

This package looks incredible. Thank you.

I want to understand something please.

You allow lidR to be extended through integrating other packages?

I want to extend the rasterize_terrain() with Shewchuk’s Triangle. It accommodates a constrained Delaunay with holes.

Would implementing it be similar to the mba example in the tutorial? Your help would be appreciated. Does that mean I would need to ‘create’ the Triangle algorithm everytime I want to execute the option?

[Jean-Romain]

According to the documentation the function Rtriangle::triangulate() Triangulate a planar straight line graph using the Triangle library. This is not what we want to generate a DTM, we want to triangulate 3D points. I don't know how you are expecting to use this package to generate a DTM. Please provide a minimal example on how to triangulate points with Rtriangle

[AdrianKriger]

Your level of service is excellent. Thank you for the extremely fast response.

I do not understand what you are saying.
A terrain is a 2.5D representation. A 2D object; embedded in 3D space, where each location (x, y) is assigned to one and only one height (z). [-> here page 1]

Are your triangulating in 3D?
While this is possible; a Delaunay triangulation is in 2D and the elevation is not taken into account. The elevation is only considered when interpolating [-> here page 42]

Perhaps I should ask a different question. How is algorithm = tin() executing a triangulation, what type of TIN is this and how are you interpolating a surface (linear or natural neighbor or etc. )?

[Jean-Romain]

tin() triangulates the 2D points and performs a linear interpolation using Z. Ok my wording what not accurate but my point remains: we need to triangulates points not straight line graph. Based on the doc of Rtriangulate I can't see any function that takes into input a set of unordered 2D coordinates. It expects an ordered set of segments representing a polygon as far as I can tell.

[AdrianKriger]

You are very kind to respond so promptly.

install.packages("devtools")
devtools::install_github("davidcsterratt/RTriangle", subdir="pkg")
library(RTriangle)

p <- pslg(P=rbind(c(0, 0), c(0, 1), c(0.5, 0.5), c(1, 1), c(1, 0)),
          S=rbind(c(1, 2), c(2, 3), c(3, 4), c(4, 5), c(5, 1)))

#Plot it
plot(p)

## Triangulate it
tp <- triangulate(p)
tp

plot(tp)

where the pslg(P, PB = NA, PA = NA, S = NA, SB = NA, H = NA) is:

P: A 2-column matrix of x-y co-ordinates of vertices. There is one row per vertex.
S: A 2-column matrix of segments in which each row is a segment. Segments are edges whose endpoints are vertices in the PSLG, and whose presence in any mesh generated from the PSLG is enforced. Each segment refers to the indices in P as the endpoints of the segment.

The segments are not necessary.

p <- pslg(P=rbind(c(0, 0), c(0, 1), c(0.5, 0.5), c(1, 1), c(1, 0)))
## Triangulate it
tp <- triangulate(p)
plot(tp)

The triangulation fed to the lidR interpolation function will be different.
Can you accommodate the library(RTriangle) triangulation?

[Jean-Romain]

Ok so we can triangulate 2D points. But you do not have any interpolation method. So we will use some lidR's internal functions for linear interpolation. However, you must be careful here, you are using my internal function with your external tools. I don't know to which extent they fit together. For exemple the triangulate function mentions it can add Steiner points. If you add points that are not a in the point cloud then it may crash, return wrong results or whatever. See the above code as a proof of concept not a code ready for production

triangulation<- function(a = NULL, q = NULL, Y = FALSE, j = FALSE, D = FALSE, S = 10000, extrapolate = knnidw(3,1,50))
{
  f = function(gnd, where)
  {
    # pslg needs a 2 coordinates object
    M = lidR:::coordinates(gnd)

    # Triangulation
    pslg = RTriangle::pslg(P = M)
    tri = RTriangle::triangulate(pslg, a, q, Y, j, D, S)

    # lidR's internal linear interpolation of triangles
    D = tri[["T"]]
    P = as.matrix(lidR:::coordinates3D(gnd))
    X = as.matrix(where)
    z = lidR:::tInterpolate(D,P,X)

    # Extrapolate beyond the convex hull
    isna <- is.na(z)
    nnas <- sum(isna)
    if (nnas > 0)
    {
      lidR.context <- "spatial_interpolation"
      where2 <- data.frame(X = where$X[isna],  Y = where$Y[isna])
      zknn <- extrapolate(gnd, where2)
      z[isna] <- zknn
    }

    return(z)
  }

  # f is a function but we can set compatible classes. Here it is an
  # algorithm for DTM
  f <- plugin_dtm(f)
  return(f)
}


library(lidR)
LASfile <- system.file("extdata", "Topography.laz", package="lidR")
las = readLAS(LASfile)

dtm = rasterize_terrain(las, 1, triangulation())
plot(dtm)
plot_dtm3d(dtm)

# generate incorrect result
dtm = rasterize_terrain(las, 1, triangulation(a = 100, S = 2))
plot(dtm)

# force 0  Steiner points (which should fall to regular Delaunay I guess)
dtm = rasterize_terrain(las, 1, triangulation(a = 100, S = 0))
plot(dtm)

[AdrianKriger]

You are very considerate. Thank you.

I need to know whether to set aside time to look into this further.

I want to exclude water (lakes, etc.) and buildings {essentially polygons} from the triangulation. Shewchucks' Triangle is extremely effective with executing this option -> [pslg(P, S, H)].

I can see:

     # pslg needs a 2 coordinates object
    M = lidR:::coordinates(gnd)

    # Triangulation
    pslg = RTriangle::pslg(P = M)

pslg takes the ground classified points but I cannot see how to pass the segments (and holes) to the pslg object. How do I pass other parameters to the pslg object please?

[Jean-Romain]

Add parameters to triangulation

triangulation<- function(PB = NA, PA = NA, S = NA, SB = NA, H = NA, a = NULL, q = NULL, Y = FALSE, j = FALSE, D = FALSE, S = 10000, extrapolate = knnidw(3,1,50))
{
  f = function(gnd, where)
  {
    # pslg needs a 2 coordinates object
    M = lidR:::coordinates(gnd)

    # Triangulation
    pslg = RTriangle::pslg(P = M, PB , PA, S, SB, H)
   [...]

[AdrianKriger]

Okay.
I will only be able to get to this in October but I'm confident the investment will bear fruit.
You have been extremely generous and helpful. Thank you.

[AdrianKriger]

It looks like the sf library already has a constrained Delaunay:

st_triangulate_constrained(po)

[Jean-Romain]

probably not one that supports hundreds of millions of points but yes you can convert laser point cloud to sf if needed with st_as_sf