Zuweisen von RGB-Werten aus dem Geotiff-Bild zu LiDAR-Daten mit R.

10

Ich habe ein Geotiff-Bild und die entsprechenden Lidar-Daten (x, y, z) in UTM-Koordinaten angegeben. Ich muss die Lidar-Daten mit den RGB-Werten aus dem Bild zusammenführen.

Das bedeutet, dass ich am Ende jeden Punkt der LiDAR-Wolkenfarbe (3D) zeichnen muss, der mit dem entsprechenden RGB-Wert aus dem Geotiff-Bild codiert ist.

Ich habe die Lidar-Daten mit QGIS in ein Shapefile konvertiert. Was soll ich als nächstes tun?

In R habe ich die plot3DFunktion ausprobiert , aber sie hat nicht funktioniert. Ich hänge das Textdokument , das Shapefile und das TIF-Bild an

Bearbeiten:

Ich habe das folgende Programm wie unten gezeigt durchgeführt:

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

Beim Versuch, das Diagramm zu zeichnen, wird jedoch der folgende Fehler angezeigt:

Fehler in [.data.frame(x @ data, i, j, ..., drop = FALSE): nicht verwendetes Argument (1)

Bearbeiten:

Ich habe das 3D-Modell ohne RGB wie unten gezeigt erhalten:

Geben Sie hier die Bildbeschreibung ein

Bibinwilson
quelle
Formdatei
bibinwilson
1
Sie verwirren Begriffe auf eine Weise, die die Frage und Ihren Code unsinnig macht. Polygone stellen diskrete Bereiche dar, während Punkte explizite x, y-Positionen sind. Es sieht so aus, als würden Sie eine Punkt-Feature-Class und kein Polygon lesen. Wenn dies der Fall ist, möchten Sie nicht "fun = mean" in der Extraktionsfunktion. Ich möchte auch darauf hinweisen, dass R nicht die ideale Software für 3D-Diagramme großer Punktwolken ist. Darüber hinaus ist Ihre Absicht für die Visualisierung in Ordnung, aber aufgrund von Parallaxenproblemen von 2D, die auf 3D-Daten projiziert werden, können Sie diese nicht analytisch verwenden.
Jeffrey Evans
Gibt es eine Möglichkeit, das Shapefile und die TIFF-Dateien zusammenzuführen, sodass ich einige andere Softwaretools verwenden kann, um sie zu zeichnen?
Bibinwilson
qustion ist einfach. Ich benötige ein 3D-Diagramm von einem RGB GEOTIFF IMAGE + XYZ-Wert.
Bibinwilson
2
Wenn Sie R nicht verwenden müssen, können Sie den PDAL-Farbfilter verwenden: pdal.io/stages/filters.colorization.html
Pete Gadomski

Antworten:

11

Vielen Dank für die Klärung Ihrer Frage, da diese bisher recht unklar war. Sie können ein Multiband-Raster mithilfe der Stack- oder Brick-Funktion im Rasterpaket lesen und die zugehörigen RGB-Werte mithilfe des Extrakts, ebenfalls aus dem Raster, einem sp SpatialPointsDataFrame-Objekt zuweisen. Das Erzwingen des data.frame-Objekts (das aus read.csv resultiert) zu einem sp-Punkt-Objekt, das zum Extrahieren übergeben werden kann, wird mit dem sp-Paket erreicht.

Das 3D-Diagramm stammt aus dem rgl-Paket. Da der Plot interaktiv ist und nicht an eine Datei übergeben wird, können Sie mit rgl.snapshot eine Datei erstellen. Die Basis-RGB-Funktion nimmt drei RGB-Werte an und erzeugt eine entsprechende einwertige R-Farbe. Durch Erstellen eines Vektors, der den Daten entspricht, können Sie ein Diagramm mit dem Argument col einfärben, ohne Farbe als tatsächliche Dimension zu definieren (was Ihre anfängliche Verwirrung zu sein schien).

Hier ist ein kurzes Dummy-Beispiel.

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

Und hier ist ein Beispiel mit den von Ihnen angegebenen Daten.

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
quelle
Ich habe den obigen Code mit den vom Poster bereitgestellten Beispieldaten ausprobiert. Es funktioniert, aber die RGB-Farben sind etwas chaotisch. Ich habe einige Dächer, die wie Straßen und umgekehrt gefärbt sind. Liegt dies wahrscheinlich an einer zu geringen Genauigkeit der Ziffern der Probe txt lidardata?
umbe1987
3

Eine Alternative zum Rendern von LiDAR-Daten und RGB-Werten in 3D ist FugroViewer .

Unten finden Sie ein Beispiel mit Beispieldaten, die sie bereitstellen. Ich habe die Datei mit Bmore_XYZIRGB.xyzdem Titel verwendet, die so aussieht:

Geben Sie hier die Bildbeschreibung ein

Wählen Sie beim Öffnen in Fugro Viewer die entsprechenden Felder aus, die in der Datei verfügbar sind (in diesem Fall eine .xyz-Datei):

Geben Sie hier die Bildbeschreibung ein

Färben Sie dann die Punkte mit den RGB-Daten aus und wählen Sie das Werkzeug aus Color Points by Encoding RGB Image Values(siehe den roten Pfeil auf dem Screenshot unten). Aktivieren Sie die 3DSchaltfläche für die 3D-Visualisierung.

Geben Sie hier die Bildbeschreibung ein

Andre Silva
quelle
3

Bearbeiten: Wie von Mathiaskopo erwähnt, verwenden neuere Versionen von LAStools Lascolor ( README ).

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

Eine andere Möglichkeit wäre, las2las wie folgt zu verwenden:

las2las -i input.las --color-source RGB_photo.tif -o output.las --file-format 1.2 --point-format 3 -v    
dmci
quelle
Die neueste Version wird mit lascolor: lascolor -i LiDAR.las -image image.tif -odix _rgb -olas
Mathiaskopo
2

Dieser Code verwendet gdal, numpy und matplotlib, um die x-, y- und z-Werte aus einem Raster zu extrahieren und ein 3D-Modell davon zu erstellen.

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

Ich habe den obigen Code mit einem Raster mit Steigungslänge (GTiff, 50 Zeilen x 50 Spalten) verwendet und das folgende Ergebnis erhalten:

Geben Sie hier die Bildbeschreibung ein

xunilk
quelle
1
Eigentlich bekomme ich das 3D-Modell. Aber ich brauche das entsprechende RGB für jedes Pixel, ich muss das aus dem GEOTiff-Bild extrahieren und das 3D-Modell
einfügen
War mein Code nützlich, um Ihr 3D-Modell zu erhalten?
Xunilk