Von XYZ in TIF mit R

Konvertierung und Verschmelzung von XYZ-Punktwolken in Raster mit R.

Die Ausgangssituation für diesen Beitrag ist das Vorliegen mehrerer Punktwolken, hier im XYZ-Format, die für Visualisierungszwecke und weitere Verarbeitung in einer gemeinsamen Rasterdatei verschmolzen werden sollen. Neben vielen anderen Möglichkeiten werden hier Punktwolken frei zum Download angeboten. In der Regel wird ein Gebiet durch mehrere Kacheln abgedeckt (siehe Kachelübersichten). Je nach Anzahl verwendeter Kacheln kann die Handhabung dieser Dateien schnell unübersichtlich werden. Zudem ist der Ressourcenbedarf von Einzelkacheln relativ hoch. Ziel ist deshalb die Umwandlung der XYZ-Dateien in Rasterdaten und die anschließende Verschmelzung zu einem Rasterdatensatz. Die Verarbeitungsprozedur wird mit Hilfe eines R-Skripts erklärt.

Zunächst wird mit der Funktion setwd() das Verzeichnis festgelegt, in dem sich die XYZ-Dateien befinden. Der Variable xyz_files werden alle Dateinamen in diesem Verzeichnis zugewiesen, deshalb sollten dort neben den XYZ-Dateien keine anderen Dateien liegen.

# read xyz point cloud
setwd("path/to/xyz_files_directory")
xyz_files = list.files(getwd())

Für die Verarbeitung wird die Bibliothek raster benötigt. Es wird über den Aufruf install.packages(„raster“) installiert. Alle Funktionalitäten werden mittels require(raster) aktiviert.

# install.packages("raster")
require(raster)

Im nächsten Schritt werden die XYZ-Dateien nacheinander mit der Funktion read.table() geladen. Dieser Funktionsaufruf muss auf die vorliegende Dateistruktur abgestimmt werden (z. B. header=TRUE wenn die erste Zeile der XYZ-Dateien eine Spaltenüberschrift enthält). Im Skript wird davon ausgegangen, dass die ersten drei Spalten jeweils Rechtswert, Hochwert und Höhenwert entsprechen (keeps <- c(„V1“, „V2“, „V3“)). Mit der Funktion rasterFromXYZ() wird anschließend jede XYZ-Datei in eine TIF-Datei umgewandelt und im oben ausgewiesenen Arbeitsverzeichnis (setwd()) gespeichert.

# convert xyz to raster (tif)
for (i in 1:length(xyz_files)) {
  print(i)
  print(length(xyz_files))
  
  xyz_file = read.table(xyz_files[i], header=F, dec=".")
  #print(head(xyz_file))
  print("loaded...")
  
  keeps <- c("V1", "V2", "V3")
  xyz_file = xyz_file[keeps]

  ras = rasterFromXYZ(xyz_file, res=c(1,1), crs=NA, digits=5)
  writeRaster(ras, paste(as.character(i), ".tif", sep=""))
  print("saved...")
}

Alle gespeicherten Dateinamen der TIF-Dateien werden im folgenden Abschnitt der Variable tif_files zugewiesen. Die Liste to_be_merged enthält nach Ausführung der for-Schleife alle TIF-Dateien und wird für die Verschmelzung als Eingangsparameter verwendet. 

# merge tifs
tif_files = list.files(getwd(), pattern=".tif$")
to_be_merged = c()
for (tif in tif_files) {
  single_raster = assign(unlist(strsplit(tif, "[.]"))[1], raster(tif))
  to_be_merged <- c(to_be_merged, single_raster)
  print(tif)
}

names(to_be_merged)[1:2] <- c('x', 'y')
to_be_merged$fun <- mean
to_be_merged$na.rm <- TRUE
merged_raster <- do.call(mosaic, to_be_merged)

Im letzten Schritt wird das verschmolzene Raster als TIF-Datei gespeichert.

# write merged raster
writeRaster(merged_raster, "merged_raster.tif")

In Folge der Verschmelzung entsteht ein großer Rasterdatensatz, dessen Ladezeit bei der Darstellung im GIS etwas länger dauern kann. Durch die zusätzliche Berechnung von Pyramiden wird die Darstellung großer Raster performanter. Das gesamte R-Skript befindet sich hier.

30. April 2020Geodaten