Adding 3d features is hard because of axis choice

[Jean-Romain]

Hi,

I was trying to render a LiDAR point cloud on top of the DTM rendered by rayshader. Using the lidR package it was expected to be easy but in practice the axes chosen in rayshader make this task particularly hard. For example the elevation is rendered on the Y axis instead of the Z axis.

Let start by making a DTM with lidR

library((lidR)
LASfile = system.file("extdata", "Topography.laz", package="lidR")
las = readLAS(LASfile)
dtm = grid_terrain(las, algorithm = tin()) # RasterLayer

Then we can compute the heightmap and calculate the surface color map.

library(rayshader)
elmat = raster_to_matrix(dtm)
elmat = t(apply(elmat, 1, rev)) # Invert of the DTM in to get back the correct topography 

x = elmat %>%
  sphere_shade(texture = "desert") %>%
  add_water(detect_water(elmat), color = "desert") %>%
  add_shadow(ray_shade(elmat, zscale = 3), 0.5) 

Now we are plotting

library(rgl)
x %>% plot_3d(elmat, zscale = 1, fov = 0, theta = 135, zoom = 0.75, phi = 45, windowsize = c(1000, 800))
axes3d(c("x-", "y-", "z-"), col = "black")
title3d(xlab = "X", ylab = "Y", zlab = "Z", col = "black")

Fisrt we can notice that the scene is rendered on (0,0) which makes a lot of sense. But then we can see that:

1. The X,Y coordinates are rendered on the X,Z axis and the elevation is rendered on the Y axis.
2. The height map has been transposed or flipped. It does not correspond to the truth. It is easy to compare with the point cloud or with a map to see that lake are not at the good place.

Then we can try to add the original point cloud. The coordinates of the points being given with a coordinate reference system they don't match with the rayshader rendering centered on (0,0) but this is not a problem because lidR already hold this case and shift coordinates for a better rendering (a story of floating point accuracy). So here we can compute the offsets to center the points on (0,0)

bbox <- extent(dtm)
xoffset = (bbox@xmax + bbox@xmin) / 2
yoffset = (bbox@ymax + bbox@ymin) / 2
offsets = c(xoffset, yoffset)

Then if the rayshader scene was rendered on X,Y,Z I could do (from lidR v3.0.0)

plot(las, add = offsets)

But it does not work because of the axes choices. So we can switch the axis manually

points3d(las$X - xoffset, las$Z, las$Y - yoffset, col = lidR:::set.colors(las$Z,  height.colors(50)))

Surprisingly it works but the scene is still inverted. What is surprising actually is the fact that points are inverted as well as if the axes themselves were inverted.

As a conclusion it is very hard to render a real scene and add custom objects. I'd like to understand your axes choices and the rendering inversion and I'd like to know if there is way to tune that easily to make this kind of rendering easier.

Thanks

ps: I know that the question is even harder with a dtm resolution different than 1 but let keep it simple with a resolution of 1.

[tylermorganwall]

The choice of the y-axis being up is a computer graphics convention, and it's the standard used in most OpenGL-based software. Using the z-axis as up is a perfectly fine choice, but it's just a different convention. As you found out, it's not too difficult to translate between the two.

The transposed geometry is due to the different data-ordering conventions used by R matrices and raster data—you either have to flip and transpose yourself, or you can use the render_points() function in rayshader, which will do it for you:

render_points(extent= bbox, long = las$X, lat = las$Y, altitude = las$Z,size=1)

[tylermorganwall]

I believe this is resolved (use the render_points() function when adding outside data, which will account for the internal orientation choices).