Przypisywanie wartości RGB z obrazu Geotiff do danych LiDAR za pomocą R.

Podałem obraz Geotiff i odpowiadające mu dane Lidara (x, y, z) we współrzędnych UTM. Muszę scalić dane Lidar z wartościami RGB z obrazu.

Oznacza to, że na koniec muszę wydrukować (3D) każdy punkt chmury LiDAR oznaczony odpowiednią wartością RGB z obrazu Geotiff.

Przekształciłem dane Lidara w plik kształtu za pomocą QGIS. Co mam teraz zrobić?

W R wypróbowałem tę plot3Dfunkcję, ale nie działała. Załączam dokument tekstowy , plik kształtu i obraz tif

Edytować:

Zrobiłem następujący program, jak pokazano poniżej:

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

Ale podczas próby wykreślenia wykresu pojawia się następujący błąd:

Błąd w [.data.frame(x @ data, i, j, ..., drop = FALSE): nieużywany argument (1)

Edytować:

Mam model 3D bez RGB, jak pokazano poniżej:

Dziękujemy za wyjaśnienie pytania, ponieważ było ono wcześniej niejasne. Możesz odczytać raster wielopasmowy za pomocą funkcji stosu lub cegły w pakiecie rastrowym i przypisać powiązane wartości RGB do obiektu SpatialPointsDataFrame za pomocą ekstrakcji, również z rastra. Przymus obiektu data.frame (wynikający z read.csv) na obiekt punktu sp, który można przekazać do wypakowania, osiąga się za pomocą pakietu sp.

Fabuła 3D pochodzi z pakietu rgl. Ponieważ wykres jest interaktywny i nie jest przekazywany do pliku, możesz utworzyć plik za pomocą rgl.snapshot. Podstawowa funkcja rgb przyjmuje trzy wartości RGB i tworzy odpowiadający jednokrotny kolor R. Tworząc wektor odpowiadający danym, możesz pokolorować wykres za pomocą argumentu col bez definiowania koloru jako rzeczywistego wymiaru (co wydawało się być początkowym zamieszaniem).

Oto szybki przykład manekina.

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

Oto działający przykład z podanymi danymi.

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

Alternatywą do renderowania danych LiDAR i wartości RGB w 3D jest FugroViewer .

Poniżej znajduje się przykład z podanymi przez nich przykładowymi danymi. Użyłem pliku zatytułowanego, Bmore_XYZIRGB.xyzktóry wygląda następująco:

Podczas otwierania w Fugro Viewer wybierz odpowiednie pola dostępne w pliku (w tym przypadku plik .xyz):

Następnie pokoloruj punkty za pomocą danych RGB, wybierając narzędzie Color Points by Encoding RGB Image Values(patrz czerwona strzałka na zrzucie ekranu poniżej). Włącz 3Dprzycisk do wizualizacji 3D.

Edycja: jak wspomniała Mathiaskopo, nowsze wersje LAStools używają lascolor ( README ).

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

Inną opcją byłoby użycie las2las w następujący sposób:

las2las -i input.las --color-source RGB_photo.tif -o output.las --file-format 1.2 --point-format 3 -v    

Ten kod używa gdal, numpy i matplotlib do wyodrębnienia wartości x, y, z rastra i uzyskania jego modelu 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

Użyłem powyższego kodu z rastrem o nachyleniu (GTiff, 50 wierszy x 50 kolumn) i uzyskałem następujący wynik: