10

I have given a Geotiff image and its corresponding Lidar data (x,y,z) in UTM coordinates. I need to merge the Lidar data with the RGB values from the image.

That means, at the end, I need to plot (3D) each point of the LiDAR cloud color coded with its corresponding RGB value from the Geotiff image.

I converted the Lidar data into a shapefile using QGIS. What should I do next?

In R, I tried the plot3D function, but, it did not work. I'm attaching the text doc, shapefile, and tif image

Edit:

I've done the following program as shown below:

require(raster) 
require(maptools)  # to take shape files
#require(car) # for scatter3D 
require(plot3Drgl)

##setwd("C:\\Users\\Bibin Wilson\\Documents\\R")
##source('Lidar.r')

data = read.csv("C:\\Users\\Bibin Wilson\\Desktop\\Lidar\\lidardata.csv")
#nr = nrow(data)
nc = ncol(data)

nr = 500

require(rgdal)
X = readGDAL("C:\\Users\\Bibin Wilson\\Desktop\\Lidar\\image.tif")

topx = 4.968622208855732e+05;
topy = 5.419739403811632e+06;

final = matrix(nrow = nr, ncol = nc+2)

for(i in 1:nr) {
 x = data[i,1]
 y = data[i,2]
 rr = round((topy-y)/0.0833)
 cc = abs(round((x-topx)/0.0833))
 if(rr == 0) {
  rr = 1
 }
 if(cc == 0) {
  cc = 1
 }
 final[i,1] = x
 final[i,2] = y
 final[i,3] = data[i,3]
 final[i,4] = rr
 final[i,5] = cc
}

for(i in 1:nr) {
 x = final[i,1]
 y = final[i,2]
 z = final[i,3]     
 rr = final[i,4]
 cc = final[i,5]
 if(rr <= 5086 && cc<=3265) {
  r = X[rr,cc,1]/255
  g = X[rr,cc,2]/255
  b = X[rr,cc,3]/255
  c = cbind(r,g,b)
  scatter3D(x,y,z,2,c)
 }
}

But while trying to plot the graph, it shows the following error:

Error in [.data.frame(x@data, i, j, ..., drop = FALSE) : unused argument (1)

Edit:

I got the 3D model without the RGB as shown below:

enter image description here

  • 1
    You are confusing terms in a way that is making the question, and your code, nonsensical. Polygons represent discrete areas whereas points are explicit x,y locations. It looks like you are reading a point feature class and not polygon. If this is the case, you do not want "fun=mean" in the extract function. I would also point out that R is not the ideal software for 3D plots of large point clouds. Additional, your intent is fine for visualization but due to parallax issues of 2D projected onto 3D data, you cannot use this analytically. – Jeffrey Evans Mar 12 '15 at 16:12
  • Is there any way to merge the shapefile and the TIFF files, so that I can use some other software tools to plot them out. – bibinwilson Mar 12 '15 at 17:08
  • qustion is simple. I need a 3D plot from one RGB GEOTIFF IMAGE + XYZ values. – bibinwilson Mar 12 '15 at 17:22
  • 2
    If you don't have to use R, you could use PDAL's colorization filter: pdal.io/stages/filters.colorization.html – Pete Gadomski Mar 16 '15 at 13:04
11

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")
  • I tried the code above with the sample data provided by the poster. It does work, but the RGB colors are a bit messy. I have som roofs coloured like streets and viceversa. Is this probabily due to too few precision in the digits of the sample txt lidardata? – umbe1987 Jan 11 '16 at 16:42
3

An alternative to render LiDAR data and RGB values in 3D is FugroViewer.

Below, there is an example with sample data they provide. I used the file entitled Bmore_XYZIRGB.xyz which looks like this:

enter image description here

When opening in Fugro Viewer select the corresponding fields available within the file (in this case, a .xyz file):

enter image description here

Then, color the points using the RGB data, selecting the tool Color Points by Encoding RGB Image Values (see the red arrow on the screenshot below). Turn on the 3D button for 3D visualization.

enter image description here

3

Edit: as mentioned by Mathiaskopo, newer versions of LAStools use lascolor (README).

lascolor -i LiDAR.las -image image.tif -odix _rgb -olas

Another option would be to use las2las as follows:

las2las -i input.las --color-source RGB_photo.tif -o output.las --file-format 1.2 --point-format 3 -v    
  • The newest version is using lascolor: lascolor -i LiDAR.las -image image.tif -odix _rgb -olas – Mathiaskopo Nov 19 '16 at 13:50
2

This code uses gdal, numpy and matplotlib for extracting the x, y, z values from a raster and to have a 3D model of it.

#!/usr/bin/env python
# -*- coding: utf-8

#Libraries
from osgeo import gdal
from os import system
import struct
import time

import numpy as np
from matplotlib.mlab import griddata
from mpl_toolkits.mplot3d.axes3d import *
from matplotlib import cm
import matplotlib.pyplot as plt

#Function to extract x,y,z values
def getCoorXYZ(band):

    # fmttypes: Byte, UInt16, Int16, UInt32, Int32, Float32 y Float64
    fmttypes = {'Byte':'B', 'UInt16':'H', 'Int16':'h', 'UInt32':'I', 'Int32':'i', 'Float32':'f', 'Float64':'d'}

    print "rows = %d columns = %d" % (band.YSize, band.XSize)

    BandType = gdal.GetDataTypeName(band.DataType)

    print "Data type = ", BandType

    x = []
    y_ = []
    z = []

    inc_x = 0

    for y in range(band.YSize):

        scanline = band.ReadRaster(0, y, band.XSize, 1, band.XSize, 1, band.DataType)
        values = struct.unpack(fmttypes[BandType] * band.XSize, scanline)

        for value in values:
            z.append(value)
            inc_x += 1
            y_.append(inc_x)
            x.append(y+1)           

        inc_x = 0

    return x, y_, z

#Program start here!

system("clear")

nameraster = str(raw_input("raster name = ? "))

start = time.time()

dataset = gdal.Open(nameraster)
band = dataset.GetRasterBand(1)

print "Processing %s" % nameraster

x,y,z = getCoorXYZ(band)

# grid 2D construction
xi = np.linspace(min(x), max(x))
yi = np.linspace(min(y), max(y))
X, Y = np.meshgrid(xi, yi)

# interpolation
Z = griddata(x, y, z, xi, yi)

#Visualization with Matplotlib
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet,linewidth=1, antialiased=True)
plt.plot

end = time.time()

time_tot = end - start

print "Total time = %.4f s" % time_tot     

plt.show() #necessary for having a static window

I used the above code with a slope lenght raster (GTiff, 50 rows x 50 columns) and I obtained the following result:

enter image description here

  • 1
    actually I'm getting the 3D model. But I need to have the corresponding RGB for each pixel, I need to extract that from the GEOTiff image and need to put in the 3D model – bibinwilson Mar 18 '15 at 12:52
  • Was my code useful to get your 3D model? – xunilk Mar 18 '15 at 13:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.