Skip to content

Instantly share code, notes, and snippets.

@wiesehahn
Created September 10, 2023 12:14
Show Gist options
  • Save wiesehahn/3f69f2b71af739169103e50ef7952f0c to your computer and use it in GitHub Desktop.
Save wiesehahn/3f69f2b71af739169103e50ef7952f0c to your computer and use it in GitHub Desktop.
Example how to render pointclouds in R with rayrender
##___________________________________________________
##
## Script name: render_pointcloud.R
##
## Purpose of script:
## render forest plot pointcloud with raytracing
##
## Author: Jens Wiesehahn
## Copyright (c) Jens Wiesehahn, 2023
## Email: wiesehahn.jens@gmail.com
##
## Date Created: 2023-09-10
##
## Notes:
## quite computationally intensive on large/dense datasets
##
##___________________________________________________
library(lidR)
library(rayshader)
library(here)
f <- system.file("extdata", "Megaplot.laz", package="lidR")
las <- readLAS(f)
nlas <- normalize_height(las, knnidw())
noground <- filter_poi(nlas, Classification != 2, Z >0.1)
noground <- unnormalize_height(noground)
# create DTM
dtm <- rasterize_terrain(las, 0.25, tin(), pkg = "raster")
# init 3D plot with DTM
elmat <- raster_to_matrix(dtm)
elmat %>%
sphere_shade(zscale = 0.25,
texture=create_texture("white","green","white","white","white")
) %>%
add_shadow(ray_shade(elmat), 0.6) %>%
plot_3d(elmat, zscale = 0.25, fov = 45, theta = 1, zoom = 0.6, phi = 89, windowsize = c(500, 500), solid=F)
# add vegetation points
render_points(extent = attr(dtm,"extent"),
lat = noground@data$Y,
long = noground@data$X,
altitude = noground@data$Z,
zscale=0.25,
offset = 0,
color=viridis::plasma(length(noground@data$Z))[rank(noground@data$Z)],
size = 3)
# get environment light scene
download.file("https://dl.polyhaven.org/file/ph-assets/HDRIs/hdr/2k/kiara_1_dawn_2k.hdr", here("kiara_1_dawn_2k.hdr"))
# render scene
render_highquality(point_radius = 3,
parallel = TRUE,
light = FALSE,
environment_light = here("kiara_1_dawn_2k.hdr"),
clear = TRUE,
width = 1000,
height = 1000,
filename = here("raytraced_pointcloud.png")
)
@hambrecht
Copy link

hambrecht commented Sep 12, 2023

Hey Jens,
Much appreciate you sharing this code example!
However, when following your example, I run into the following issue when executing render_highquality:

Error in rayvertex::validate_mesh(mesh) : 
  !any(is.na(normals)) is not TRUE

I get the same error when testing it with my data.
Any ideas?

UPDATE:
My current suspicion is that it is caused by NA values in the DTM. Potential fix:

library(imputeTS)
dtm_interpolated <- na.interpolation(dtm, option = "linear")

@wiesehahn
Copy link
Author

It works for me, so I guess its your data, otherwise you might open an issue at https://github.com/tylermorganwall/rayrender

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment