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

10

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:

wprowadź opis zdjęcia tutaj

bibinwilson
źródło
1
Mylisz terminy w sposób, który sprawia, że ​​pytanie, a Twój kod jest nonsensowny. Wieloboki reprezentują dyskretne obszary, podczas gdy punkty są wyraźnymi lokalizacjami x, y. Wygląda na to, że czytasz klasę punktów, a nie wielokąt. W takim przypadku nie chcesz „fun = mean” w funkcji wyodrębniania. Chciałbym również zauważyć, że R nie jest idealnym oprogramowaniem do wykresów 3D dużych chmur punktów. Dodatkowo, twoje zamiary są dobre do wizualizacji, ale z powodu problemów z paralaksą 2D wyświetlanych na danych 3D, nie możesz tego użyć analitycznie.
Jeffrey Evans,
Czy jest jakiś sposób scalenia pliku shapefile i plików TIFF, abym mógł użyć innych narzędzi programowych do ich wykreślenia.
bibinwilson
pytanie jest proste. Potrzebuję wykresu 3D z jednego obrazu GEOTIFF RGB + wartości XYZ.
bibinwilson
2
Jeśli nie musisz używać R, możesz użyć filtru koloryzacji PDAL
Pete Gadomski

Odpowiedzi:

11

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")
Jeffrey Evans
źródło
Próbowałem powyższego kodu z przykładowymi danymi dostarczonymi przez plakat. Działa, ale kolory RGB są nieco niechlujne. Mam dachy w kolorze jak ulice i viceversa. Czy jest to prawdopodobnie spowodowane zbyt małą precyzją cyfr próbki TXT?
umbe1987
3

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:

wprowadź opis zdjęcia tutaj

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

wprowadź opis zdjęcia tutaj

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.

wprowadź opis zdjęcia tutaj

Andre Silva
źródło
3

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    
dmci
źródło
Najnowsza wersja używa lascolor: lascolor -i LiDAR.las -image image.tif -odix _rgb -olas
Mathiaskopo
2

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:

wprowadź opis zdjęcia tutaj

Xunilk
źródło
1
właściwie dostaję model 3D. Ale muszę mieć odpowiedni RGB dla każdego piksela, muszę go wyodrębnić z obrazu GEOTiff i muszę umieścić w modelu 3D
bibinwilson
Czy mój kod był przydatny do uzyskania modelu 3D?
xunilk