Suggestion: use raster/quadmesh forms for DTM

[mdsumner]

I'm delighted to see the terrain-extraction facilities, this is really awesome! I see that lasmetrics from grid_terrain is a fully-expanded data frame, so I tried converting to raster, and then to the quad-mesh structures used by rgl.

If having this in lidR is of interest I'd be happy to help, the conversion to rgl is pretty straightforward but is a bit mysterious. This is the perfect application for using it to save on the size of the gridded terrain - and the quads can be reprojected without problems, and they can also be used to to texture map on overlay images. I don't think lidR already has this kind of example?

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

dtm = grid_terrain(lidar, method = "knnidw", k = 6)
## convert to regular raster, removing unnessary coordinate storage
library(raster)
r <- rasterFromXYZ(dtm)
pryr::object_size(dtm)  # 2.52 MB
pryr::object_size(r)    # 850 kB

## convert to quadmesh (for rgl)
library(quadmesh)
qm <- quadmesh(r)

library(rgl)
cols <- viridis::viridis(100)
shade3d(qm, col = cols[scales::rescale(qm$vb[3,qm$ib]) * (length(cols) -1) + 1])

[mdsumner]

Here's a texture mapping example, it's hard to tell if this is actually right but I think it is. I'm assuming the Lidar example data is in the same UTM 17 as the lakes shapefile is.

library(lidR)
library(raster)
library(quadmesh)
library(dismo)
library(rgl)
LASfile <- system.file("extdata", "Topography.laz", package="lidR")
lidar = readLAS(LASfile)
shapefile_dir <- system.file("extdata", package = "lidR")

dtm =   grid_terrain(lidar, res = 1, method = "knnidw", k = 6L)
lakes = rgdal::readOGR(shapefile_dir, "lake_polygons_UTM17")

r <- rasterFromXYZ(dtm, crs = projection(lakes))
qm <- quadmesh(r)

## download a google satellite image with dismo
gm <- gmap(x = lakes, type = "satellite", scale = 2)

## 1. Create PNG for texture
# we need RGB expanded (gmap gives a palette)
rgb1 <- col2rgb(gm@legend@colortable)
img <- brick(gm, gm, gm)
cells <- values(gm) + 1
img <- setValues(img, cbind(rgb1[1, cells], rgb1[2, cells], rgb1[3, cells]))
## finally, create RGB PNG image to act as a texture image
pngfile <- sprintf("%s.png", tempfile())
rgdal::writeGDAL(as(img, "SpatialGridDataFrame"), pngfile, drivername = "PNG", type = "Byte", mvFlag = 255)

## 2. Remap the image coordinates (Mercator) onto (longlat) then 
## into the local UTM projection and then convert to PNG [0, 1, 0, 1]
tcoords <- xyFromCell(setExtent(img, extent(0, 1, 0, 1)), 
                      cellFromXY(img, rgdal::project(
                                 rgdal::project(t(qm$vb[1:2, ]), projection(lakes), inv = TRUE), 
                                 projection(img))))

shade3d(qm, col = "white", texture = pngfile, texcoords = tcoords[qm$ib, ])

There's unfortunately no way to pass a texture to rgl without an external PNG file, hence all the munging to get that image in the right form.

I don't mind if you don't want the extra work and ignore this issue, but I'd like to leave this example here because lidR is what I'll be using when I explore this more deeply. Thanks!

[Jean-Romain]

Hi Michael,

lidR does not return expended data.frame but expended data.table. And this tiny difference change everything internally for efficiency, computation time and memory optimization even if a raster is lighter than an expended table. I have many reasons for returning expended table instead of raster. However I know that people usually need to convert to raster so there is a function as.raster in the package.

dtm =   grid_terrain(lidar, res = 1, method = "knnidw", k = 6L)
rdtm = as.raster(dtm)

Regarding quadmesh there is a function plot3d in the package. I assume rgl does the same on the fly. Do you really need to get a variable inside the R environment?

plot3d(dtm)
plot3d(rdtm) # should work as well

Concerning your second post, what exactly would you like to have in the package? I like the example but I think the function lasclassify does almost the same job (see wiki exemples) but at the point cloud level. Indeed lidR is more point-cloud oriented. I believe that the tools to manipulate rasters, vector layers, maps and so on already exist in many other very good packages and I not really keen to redo what already exists.

I'm open to any proposal. As you know because we already discussed about that. I'm not a geospatial guy. I'm not necessarily familiar with usual process and data structure used in geospatial processes so I'm really keen to listen to user advices. However, in my opinion there are also important points to consider such as:

• lidR is designed to be be the most efficient as possible (within the limits of my knowledge obviously), often against the R language itself and that is why many thing rely on data.table. And the natural outputs of data.table processes are data.table objects. That explain a lot of things.
• lidR was initially designed to my own work which implied unconventional manipulations. But since two years I improve it only to the benefit of the community, I do not really use it myself. Some of the function are designed for consistency with older code. I do agree that grid_terrain could return a raster but this would no be consistent with other functions. I prefer to use as.raster and it's okay for the extra memory used. Actually I'm much more concerned with useless in-memory copies of the point clouds usually made by regular users (usually without knowing it) than with this tiny problem.

[mdsumner]

Hi thanks for the detailed reply, lots of good information. I'll explore more, no definite request yet. (Having ground detection so accessible is really handy!)

[Jean-Romain]

No problem. If you have feature request do not hesitate. For PR too.