Atribuindo valores RGB da imagem do Geotiff aos dados LiDAR, usando R

10

Eu dei uma imagem Geotiff e seus dados Lidar correspondentes (x, y, z) nas coordenadas UTM. Preciso mesclar os dados do Lidar com os valores RGB da imagem.

Isso significa que, no final, eu preciso plotar (3D) cada ponto da nuvem LiDAR codificada por cores com seu valor RGB correspondente da imagem Geotiff.

Eu converti os dados do Lidar em um shapefile usando o QGIS. O que eu devo fazer a seguir?

Em R, tentei a plot3Dfunção, mas não funcionou. Estou anexando o documento doc , shapefile e tif image

Editar:

Eu fiz o seguinte programa, como mostrado abaixo:

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)
 }
}

Mas, ao tentar plotar o gráfico, ele mostra o seguinte erro:

Erro em [.data.frame(x @ dados, i, j, ..., drop = FALSE): argumento não utilizado (1)

Editar:

Eu consegui o modelo 3D sem o RGB, como mostrado abaixo:

insira a descrição da imagem aqui

bibinwilson
fonte
1
Você está confundindo termos de uma maneira que torna a pergunta e seu código sem sentido. Os polígonos representam áreas discretas, enquanto os pontos são locais x, y explícitos. Parece que você está lendo uma classe de recurso de ponto e não um polígono. Se for esse o caso, você não deseja "fun = mean" na função extrair. Eu também apontaria que R não é o software ideal para gráficos 3D de grandes nuvens de pontos. Além disso, sua intenção é boa para visualização, mas devido a problemas de paralaxe de 2D projetados em dados 3D, você não pode usá-lo analiticamente.
21415 Jeffrey Evans
Existe alguma maneira de mesclar os arquivos shapefile e TIFF, para que eu possa usar outras ferramentas de software para plotá-los.
bibinwilson
pergunta é simples. Preciso de um gráfico 3D a partir de um valor RGB GEOTIFF IMAGE + XYZ.
bibinwilson
2
Se você não tem que usar R, você poderia usar filtro de PDAL colorização: pdal.io/stages/filters.colorization.html
Pete Gadomski

Respostas:

11

Obrigado por esclarecer sua pergunta, pois antes ela não era clara. Você pode ler uma varredura multibanda usando a função stack ou brick no pacote raster e atribuir os valores RGB associados a um objeto SpatialPointsDataFrame sp usando extração, também da varredura. A coerção do objeto data.frame (que resulta de read.csv) para um objeto de ponto sp, que pode ser transmitido para extração, é obtida usando o pacote sp.

A plotagem 3D vem do pacote rgl. Como a plotagem é interativa e não passada para um arquivo, é possível criar um arquivo usando rgl.snapshot. A função rgb base utiliza três valores RGB e cria uma cor R de valor único correspondente. Ao criar um vetor, correspondente aos dados, é possível colorir uma plotagem usando o argumento col sem definir a cor como uma dimensão real (que parecia ser sua confusão inicial).

Aqui está um exemplo rápido e fictício.

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")

E, aqui está um exemplo trabalhado com os dados que você forneceu.

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")
Jeffrey Evans
fonte
Eu tentei o código acima com os dados de amostra fornecidos pelo pôster. Funciona, mas as cores RGB são um pouco confusas. Tenho alguns telhados coloridos como ruas e vice-versa. Provavelmente, isso se deve a pouca precisão nos dígitos da amostra txt lidardata?
umbe1987
3

Uma alternativa para renderizar dados LiDAR e valores RGB em 3D é o FugroViewer .

Abaixo, há um exemplo com dados de amostra que eles fornecem. Eu usei o arquivo intitulado Bmore_XYZIRGB.xyzque se parece com isso:

insira a descrição da imagem aqui

Ao abrir no Fugro Viewer, selecione os campos correspondentes disponíveis no arquivo (neste caso, um arquivo .xyz):

insira a descrição da imagem aqui

Em seguida, pinte os pontos usando os dados RGB, selecionando a ferramenta Color Points by Encoding RGB Image Values(veja a seta vermelha na captura de tela abaixo). Ligue o 3Dbotão para visualização 3D.

insira a descrição da imagem aqui

Andre Silva
fonte
3

Edit: como mencionado por Mathiaskopo, as versões mais recentes do LAStools usam lascolor ( README ).

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

Outra opção seria usar las2las da seguinte maneira:

las2las -i input.las --color-source RGB_photo.tif -o output.las --file-format 1.2 --point-format 3 -v    
dmci
fonte
A versão mais recente está usando lascolor: lascolor -i LiDAR.las -image image.tif -odix _rgb -olas
Mathiaskopo
2

Esse código usa gdal, numpy e matplotlib para extrair os valores x, y, z de uma varredura e ter um modelo 3D.

#!/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

Usei o código acima com uma varredura de comprimento de inclinação (GTiff, 50 linhas x 50 colunas) e obtive o seguinte resultado:

insira a descrição da imagem aqui

xunilk
fonte
1
na verdade, estou recebendo o modelo 3D. Mas eu preciso ter o RGB correspondente para cada pixel, eu preciso extrato que a partir da imagem GeoTIFF e necessidade de colocar no modelo 3D
bibinwilson
Meu código foi útil para obter seu modelo 3D?
xunilk