Asignación de valores RGB de imagen Geotiff a datos LiDAR, usando R

10

He dado una imagen de Geotiff y sus datos Lidar correspondientes (x, y, z) en coordenadas UTM. Necesito fusionar los datos Lidar con los valores RGB de la imagen.

Eso significa que, al final, necesito trazar (3D) cada punto del color de la nube LiDAR codificado con su valor RGB correspondiente de la imagen Geotiff.

Convertí los datos de Lidar en un shapefile usando QGIS. ¿Qué debería hacer después?

En R, probé la plot3Dfunción, pero no funcionó. Adjunto el documento de texto , el archivo de forma y la imagen tif

Editar:

Hice el siguiente programa como se muestra a continuación:

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

Pero al intentar trazar el gráfico, muestra el siguiente error:

Error en [.data.frame(x @ data, i, j, ..., drop = FALSE): argumento no utilizado (1)

Editar:

Obtuve el modelo 3D sin RGB como se muestra a continuación:

ingrese la descripción de la imagen aquí

bibinwilson
fuente
1
Estás confundiendo términos de una manera que hace que la pregunta y tu código no tengan sentido. Los polígonos representan áreas discretas, mientras que los puntos son ubicaciones explícitas x, y. Parece que está leyendo una clase de entidad de puntos y no un polígono. Si este es el caso, no desea "diversión = media" en la función de extracción. También quisiera señalar que R no es el software ideal para gráficos 3D de nubes de puntos grandes. Además, su intención está bien para la visualización, pero debido a problemas de paralaje de 2D proyectado en datos 3D, no puede usar esto analíticamente.
Jeffrey Evans
¿Hay alguna manera de fusionar el archivo de forma y los archivos TIFF, de modo que pueda usar otras herramientas de software para trazarlos?
bibinwilson
La pregunta es simple. Necesito un diagrama 3D de uno de los valores RGB GEOTIFF IMAGE + XYZ.
bibinwilson
2
Si no tiene que usar R, puede usar el filtro de colorización de PDAL
Pete Gadomski

Respuestas:

11

Gracias por aclarar su pregunta, ya que anteriormente no estaba clara. Puede leer un ráster multibanda usando la función de pila o ladrillo en el paquete ráster y asignar los valores RGB asociados a un objeto sp SpatialPointsDataFrame usando extracto, también desde el ráster. La coerción del objeto data.frame (que resulta de read.csv) a un objeto de punto sp, que se puede pasar a extraer, se logra utilizando el paquete sp.

La trama 3D proviene del paquete rgl. Dado que el gráfico es interactivo y no se pasa a un archivo, puede crear un archivo usando rgl.snapshot. La función base rgb toma tres valores RGB y crea un color R de valor único correspondiente. Al crear un vector, correspondiente a los datos, puede colorear un gráfico utilizando el argumento col sin definir el color como una dimensión real (que parecía ser su confusión inicial).

Aquí hay un ejemplo ficticio rápido.

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

Y, aquí hay un ejemplo trabajado con los datos que proporcionó.

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
fuente
Probé el código anterior con los datos de muestra proporcionados por el póster. Funciona, pero los colores RGB son un poco desordenados. Tengo algunos techos coloreados como calles y viceversa. ¿Es esta probabilidad debido a muy poca precisión en los dígitos de la muestra tidar lidardata?
umbe1987
3

Una alternativa para representar datos LiDAR y valores RGB en 3D es FugroViewer .

A continuación, hay un ejemplo con datos de muestra que proporcionan. Usé el archivo titulado Bmore_XYZIRGB.xyzque se ve así:

ingrese la descripción de la imagen aquí

Al abrir en Fugro Viewer, seleccione los campos correspondientes disponibles dentro del archivo (en este caso, un archivo .xyz):

ingrese la descripción de la imagen aquí

Luego, colorea los puntos usando los datos RGB, seleccionando la herramienta Color Points by Encoding RGB Image Values(ver la flecha roja en la captura de pantalla a continuación). Active el 3Dbotón para visualización en 3D.

ingrese la descripción de la imagen aquí

Andre Silva
fuente
3

Editar: como lo menciona Mathiaskopo, las versiones más nuevas de LAStools usan lascolor ( README ).

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

Otra opción sería usar las2las de la siguiente manera:

las2las -i input.las --color-source RGB_photo.tif -o output.las --file-format 1.2 --point-format 3 -v    
dmci
fuente
La versión más nueva está usando lascolor: lascolor -i LiDAR.las -image image.tif -odix _rgb -olas
Mathiaskopo
2

Este código usa gdal, numpy y matplotlib para extraer los valores x, y, z de un ráster y tener un modelo 3D del mismo.

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

Utilicé el código anterior con un ráster de longitud de pendiente (GTiff, 50 filas x 50 columnas) y obtuve el siguiente resultado:

ingrese la descripción de la imagen aquí

xunilk
fuente
1
De hecho, estoy obteniendo el modelo 3D. Pero necesito tener el RGB correspondiente para cada píxel, necesito extraerlo de la imagen de GEOTiff y ponerlo en el modelo 3D
bibinwilson
¿Fue útil mi código para obtener su modelo 3D?
xunilk