Affectation de valeurs RVB de l'image Geotiff aux données LiDAR, à l'aide de R

10

J'ai donné une image Geotiff et ses données Lidar correspondantes (x, y, z) en coordonnées UTM. J'ai besoin de fusionner les données Lidar avec les valeurs RVB de l'image.

Cela signifie qu'à la fin, je dois tracer (3D) chaque point de la couleur du nuage LiDAR codé avec sa valeur RVB correspondante à partir de l'image Geotiff.

J'ai converti les données Lidar en un fichier de formes à l'aide de QGIS. Que devrais-je faire ensuite?

Dans R, j'ai essayé la plot3Dfonction, mais cela n'a pas fonctionné. J'attache le document texte , le fichier de formes et l' image tif

Éditer:

J'ai fait le programme suivant comme indiqué ci-dessous:

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

Mais en essayant de tracer le graphique, il montre l'erreur suivante:

Erreur dans [.data.frame(x @ data, i, j, ..., drop = FALSE): argument inutilisé (1)

Éditer:

J'ai obtenu le modèle 3D sans RVB comme indiqué ci-dessous:

entrez la description de l'image ici

bibinwilson
la source
1
Vous confondez les termes d'une manière qui rend la question et votre code insensés. Les polygones représentent des zones discrètes tandis que les points sont des emplacements x, y explicites. Il semble que vous lisiez une classe d'entités ponctuelles et non un polygone. Si c'est le cas, vous ne voulez pas "fun = mean" dans la fonction d'extraction. Je voudrais également souligner que R n'est pas le logiciel idéal pour les tracés 3D de grands nuages ​​de points. De plus, votre intention est bonne pour la visualisation, mais en raison de problèmes de parallaxe de la projection 2D sur des données 3D, vous ne pouvez pas l'utiliser de manière analytique.
Jeffrey Evans
Existe-t-il un moyen de fusionner le fichier de formes et les fichiers TIFF, afin que je puisse utiliser d'autres outils logiciels pour les tracer.
bibinwilson
qustion est simple. J'ai besoin d'un tracé 3D à partir d'une valeur RGB GEOTIFF IMAGE + XYZ.
bibinwilson
2
Si vous n'avez pas à utiliser R, vous pouvez utiliser le filtre de colorisation de PDAL
Pete Gadomski

Réponses:

11

Merci d'avoir clarifié votre question car elle n'était pas assez claire auparavant. Vous pouvez lire un raster multibande à l'aide de la fonction pile ou brique du package raster et affecter les valeurs RVB associées à un objet Sp SpatialPointsDataFrame à l'aide de l'extraction, également à partir du raster. La contrainte de l'objet data.frame (qui résulte de read.csv) sur un objet point sp, qui peut être passé pour être extrait, est obtenue à l'aide du package sp.

L'intrigue 3D provient du package rgl. Le tracé étant interactif et non transmis à un fichier, vous pouvez créer un fichier à l'aide de rgl.snapshot. La fonction RVB de base prend trois valeurs RVB et crée une couleur R à valeur unique correspondante. En créant un vecteur, correspondant aux données, vous pouvez colorier un tracé en utilisant l'argument col sans définir la couleur comme une dimension réelle (ce qui semblait être votre confusion initiale).

Voici un exemple factice rapide.

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

Et, voici un exemple concret avec les données que vous avez fournies.

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
la source
J'ai essayé le code ci-dessus avec les exemples de données fournis par l'affiche. Cela fonctionne, mais les couleurs RVB sont un peu désordonnées. J'ai des toits colorés comme des rues et vice-versa. Est-ce probablement dû à trop peu de précision dans les chiffres de l'échantillon lidardata txt?
umbe1987
3

Une alternative pour rendre les données LiDAR et les valeurs RVB en 3D est FugroViewer .

Ci-dessous, un exemple avec des exemples de données qu'ils fournissent. J'ai utilisé le fichier intitulé Bmore_XYZIRGB.xyzqui ressemble à ceci:

entrez la description de l'image ici

Lors de l'ouverture dans Fugro Viewer, sélectionnez les champs correspondants disponibles dans le fichier (dans ce cas, un fichier .xyz):

entrez la description de l'image ici

Ensuite, coloriez les points à l'aide des données RVB, en sélectionnant l'outil Color Points by Encoding RGB Image Values(voir la flèche rouge sur la capture d'écran ci-dessous). Activez le 3Dbouton pour la visualisation 3D.

entrez la description de l'image ici

Andre Silva
la source
3

Edit: comme mentionné par Mathiaskopo, les nouvelles versions de LAStools utilisent lascolor ( README ).

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

Une autre option serait d'utiliser las2las comme suit:

las2las -i input.las --color-source RGB_photo.tif -o output.las --file-format 1.2 --point-format 3 -v    
dmci
la source
La dernière version utilise lascolor: lascolor -i LiDAR.las -image image.tif -odix _rgb -olas
Mathiaskopo
2

Ce code utilise gdal, numpy et matplotlib pour extraire les valeurs x, y, z d'un raster et pour en avoir un modèle 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

J'ai utilisé le code ci-dessus avec un raster de longueur de pente (GTiff, 50 lignes x 50 colonnes) et j'ai obtenu le résultat suivant:

entrez la description de l'image ici

xunilk
la source
1
en fait, je reçois le modèle 3D. Mais je dois avoir le RVB correspondant pour chaque pixel, je dois l'extraire de l'image GEOTiff et je dois le mettre dans le modèle 3D
bibinwilson
Mon code était-il utile pour obtenir votre modèle 3D?
xunilk