Thank you for clarifying your question as it was previously quite unclear. You can read a multiband raster using the stack or brick function in the raster package and assign the associated RGB values to an sp SpatialPointsDataFrame object using extract, also from raster. Coercion of the data.frame object (which results from read.csv) to an sp point object,that can be passed to extract, is achieved using the sp package.
The 3D plot comes from the rgl package. Since the plot is interactive and not passed to a file, you can create a file using rgl.snapshot. The base rgb function takes three RGB values and creates a corresponding single-value R color. By creating a vector, corresponding to the data, you can color a plot using the col argument without defining color as an actual dimension (which seemed to be your initial confusion).
Here is a quick dummy example.
require(rgl)
require(sp)
n=100
# Create a dummy datafame object with x,y,z values
lidar <- data.frame(x=runif(n,1,10), y=runif(n,1,10), z=runif(n,0,50))
coordinates(lidar) <- ~x+y
# Add dummy RGB values
lidar@data <- data.frame(lidar@data, red=round(runif(n,0,255),0), green=round(runif(n,0,255),0),
blue=round(runif(n,0,255),0))
# Create color vector using rgb values
cols <- rgb(lidar@data[,2:4], maxColorValue = 255)
# Interactive 3D plot
plot3d(coordinates(lidar)[,1],coordinates(lidar)[,2],lidar@data[,"z"], col=cols,
pch=18, size=0.75, type="s", xlab="x", ylab="x", zlab="elevation")
And, here is a worked example with the data you provided.
require(raster)
require(rgl)
setwd("D:/TMP")
# read flat file and assign names
lidar <- read.table("lidar.txt")
names(lidar) <- c("x","y","z")
# remove the scatter outlier(s)
lidar <- lidar[lidar$z >= 255 ,]
# Coerce to sp spatialPointsDataFrame object
coordinates(lidar) <- ~x+y
# subsample data (makes more tractable but not necessary)
n=10000
lidar <- lidar[sample(1:nrow(lidar),n),]
# Read RGB tiff file
img <- stack("image.tif")
names(img) <- c("r","g","b")
# Assign RGB values from raster to points
lidar@data <- data.frame(lidar@data, extract(img, lidar))
# Remove NA values so rgb function will not fail
na.idx <- unique(as.data.frame(which(is.na(lidar@data), arr.ind = TRUE))[,1])
lidar <- lidar[-na.idx,]
# Create color vector using rgb values
cols <- rgb(lidar@data[,2:4], maxColorValue = 255)
# Interactive 3D plot
plot3d(coordinates(lidar)[,1],coordinates(lidar)[,2],lidar@data[,"z"], col=cols,
pch=18, size=0.35, type="s", xlab="x", ylab="x", zlab="elevation")
Looks like you haven't defined a variable called image - which is what you pass along to the visualize() function. Try changing your visualize line to:
var imageRGB = dataset.median().visualize({bands: ['B4', 'B3', 'B2'], min: 0, max: 2500});
print("rgb",imageRGB);
You’ll also need to update the geometry in your export because ‘image’ doesn’t exist:
Export.image.toDrive({
scale: 30,
image: imageRGB,
description: 'example_image_RGB',
fileFormat: 'GeoTIFF',
region: geometry});
Best Answer
There is no specific GIS software for doing this: most will handle the RGB image and the Lidar data.
Basically, NDVI is (NIR - RED)/(NIR + RED). Most of the time, aerial Lidar gives you the NIR value (to be checked in metadata) and the first band of your RGB image gives you the RED value. Just make sure that your data are calibrated to reflectance (or, if not possible, normalised in the same range of value -> this will give you a good discrimination, but no absolute threshold can then be used).
As a remark, Lidar does not always provide the intensity information. Sometimes only the distance information is available.