Übung 4
Der
Normalized Burn Ratio (NBR) ist ein vegetationsspezifischer Index in der Fernerkundung, der speziell dazu entwickelt wurde Brandflächen zu identifizieren und die Schwere von Brandschäden (Burn Severity) zu bewerten. Er nutzt Satellitendaten (wie Landsat oder Sentinel-2), um Vegetationsveränderungen vor und nach einem Feuer zu vergleichen.
Ziel ist es, das Ergebnis der Klassifikation sowie die finale Karte zu verbessern. Führen Sie dazu folgende Analysen durch:
Berechnen Sie den Normalised Burn Ratio (NBR) (siehe z. B. Beetz et al. 2024) auf die gleiche Weise, wie zuvor der NDVI berechnet wurde.
Untersuchen Sie die Korrelation des NBR mit anderen Merkmalen und analysieren Sie, wie gut sich der NBR zur Trennung der Störungsklassen eignet.
Trainieren Sie das Random-Forest-Modell erneut, diesmal unter zusätzlicher Verwendung des NBR.
Vergleichen Sie die Leistungsmetriken des Modells mit den vorherigen Ergebnissen (also ohne NBR).
Beschreiben Sie, ob der Einsatz des NBR die Klassifikation verbessert hat.
Wenden Sie das trainierte Modell auf das vollständige Kompositbild an.
Aktualisieren Sie die CLC-Karte und diskutieren Sie die Ergebnisse.
Laden Sie zusätzlich die Landsat-8- und Sentinel-2-Bilder vom 10. September 2021 und erstellen Sie ein Komposit für 2021.
Berechnen Sie den Difference Normalised Burn Ratio (dNBR) aus dem NBR von 2021 und 2022 (siehe z. B. Beetz et al. 2024).
Wiederholen Sie anschließend die Analyse aus 7.1:
Untersuchen Sie die Modellgüte beim Training.
Analysieren Sie die finale Karte.
Hat die Verwendung von dNBR die Ergebnisse verbessert?
Angfangen allgemeinen Code zu schreiben, aber aufgehört und auf Aufgaben konzentriert.
library(terra)
library(RColorBrewer)
library(randomForest)
library(caret)
library(dplyr)
————————————————————————————————
l9_050922.r <- rast("data/l9_050922.tif")
plotRGB(l9_050922.r, stretch="lin", main="Landsat 9 (2022-09-05)")
plot(l9_050922.r)
crs <- crs(l9_050922.r)
crs # UTM-Zone 33N
s2a_050922.r <- rast("data/s2a050922_clouds.tif")
plotRGB(s2a_050922.r, stretch="lin", main="Sentinel-2a (2022-09-05)")
s2b_070922.r <- rast("data/s2b070922_clouds.tif")
plotRGB(s2b_070922.r, stretch="lin", main="Sentinel-2b (2022-09-07)")
— —- —- —- —- - — - —- —- — - - —- —— - —— ——- - ——- ——
# 1 Pakete installieren und laden --------------------------------------------------
# Notwendige Pakete installieren
source("install_packages.R")
# Notwendige Bibliotheken laden
library(terra) # für die Arbeit mit Rasterdaten
library(RColorBrewer) # für Farbpaletten
library(randomForest) # der Random-Forest-Klassifikator
library(caret) # für Training, Validierung und Regressions-/Klassifikationsalgorithmen
library(dplyr) # zur Arbeit mit und Manipulation von Data Frames
# 2 Daten einlesen und vorbereiten --------------------------------------------------
# 2.1 Satellitendaten/Karten einlesen --------------------------------------------
# Landsat-9 vom 5. September 2022
# Das Koordinatenreferenzsystem der Landsat-Daten wird für alle anderen Datensätze verwendet
# Sentinel-2a vom 5. September 2022
# Sentinel-2b vom 7. September 2022
# 2.2 Ein Kompositbild aus Landsat + Sentinel-2 erstellen -----------------------
# Das Landsat-Bild als Basis für das Komposit verwenden
composite_0922.r <- l9_050922.r
# Sentinel-2-Daten auf die Landsat-Daten resamplen
s2a_050922.r <- resample(s2a_050922.r, composite_0922.r)
s2b_070922.r <- resample(s2b_070922.r, composite_0922.r)
# Ein Komposit erstellen, indem der Mittelwert aller drei Datenquellen gebildet und kombiniert wird
composite_0922.r$BLUE <- mean(c(l9_050922.r$BLUE, s2a_050922.r$Blue, s2b_070922.r$Blue), na.rm=TRUE)
composite_0922.r$GREEN <- mean(c(l9_050922.r$GREEN, s2a_050922.r$Green, s2b_070922.r$Green), na.rm=TRUE)
composite_0922.r$RED <- mean(c(l9_050922.r$RED, s2a_050922.r$Red, s2b_070922.r$Red), na.rm=TRUE)
composite_0922.r$NIR <- mean(c(l9_050922.r$NIR, s2a_050922.r$NIR, s2b_070922.r$NIR), na.rm=TRUE)
composite_0922.r$SWIR1 <- mean(c(l9_050922.r$SWIR1, s2a_050922.r$SWIR1, s2b_070922.r$SWIR1), na.rm=TRUE)
composite_0922.r$SWIR2 <- mean(c(l9_050922.r$SWIR2, s2a_050922.r$SWIR2, s2b_070922.r$SWIR2), na.rm=TRUE)
plot(composite_0922.r)
plotRGB(composite_0922.r, stretch="lin", main="L9 + S2 Komposit (2022-09)")
# 2.3 Corine Land Cover+ Backbone-Karte einlesen und reprojizieren --------------
# Rasterkarte einlesen
clc.r <- rast("data/CLCplus_RASTER_2021_010m_03035/CLCplus_RASTER_2021_010m_03035.tif")
# Legende einlesen
clc.lgd <- read.csv("data/CLCplus_RASTER_2021_010m_03035/CLCplus_RASTER_2021_010m_03035_legend.csv", sep=",")
clc.lgd
# Legende dem Raster zuweisen
levels(clc.r) <- clc.lgd
coltab(clc.r) <- clc.lgd[,-2]
plot(clc.r)
# CLC-Karte projizieren
clc.r <- project(clc.r, crs, method="near")
# Auf das Ausmaß des Untersuchungsgebiets zuschneiden
clc.r <- crop(clc.r, ext(composite_0922.r))
# 2.4 Felddaten einlesen und in räumliche Objekte umwandeln -----------------------
# Felddaten zu Störungen und Composite Burn Index (CBI) einlesen
field.df <- read.csv("data/Beetz_etal_field_surveys2022-2023.csv")
field.df$disturbance <- as.factor(field.df$disturbance)
summary(field.df)
plot(field.df) # Dies plottet ein normales data.frame
# Die CBI-Felddaten von einem data.frame in einen räumlichen Vektor umwandeln
# Die Funktion erkennt automatisch die Spalten lon und lat als Koordinaten
field.v <- vect(field.df, crs="+proj=longlat +datum=WGS84")
plot(field.v, "disturbance") # Plottet nun den räumlichen Vektor in geografischen Koordinaten
# Projektion der Felddaten vom geografischen Koordinatensystem in das
# Koordinatensystem der Satellitendaten (Landsat-9)
field.v <- project(field.v, crs)
plot(field.v, "disturbance") # Plottet den räumlichen Vektor in UTM-Koordinaten
# 2.5 Daten gemeinsam in einer Karte darstellen ------------------------------------
# Felddaten über die Landbedeckungskarte legen
plot(clc.r, plg=list(title="Landbedeckung"), mar=c(2, 2, 2, 12))
plot(field.v, "disturbance", add=TRUE, col=c("black", "darkblue", "grey"),
plg=list(x="bottomright", bg="white", bty="o", title="Störung", cex=0.8))
# Landsat-Bild im Bereich der CBI-Felddaten darstellen
plotRGB(composite_0922.r, r=4, g=3, b=2, stretch="lin", axes=TRUE,
ext=ext(440000, 455000, 5632000, 5642000), mar=c(2, 2, 2, 1))
plot(field.v, "CBI", add=TRUE, col=brewer.pal(9, "Greys"),
plg=list(x="bottomleft", bg="white", bty="o", title="Brandschwere (CBI)", cex=0.8))
# Kannst du die Brandfläche des Feuers von 2022 erkennen?
# 3 Trainingsdatensatz für den Klassifikator erstellen -------------------------------
# Eine Kopie der Felddaten erstellen
data.df <- field.df
# Die Bänder des Landsat-Bildes zum Felddatensatz hinzufügen
xy <- crds(field.v) # Koordinaten aus den Felddaten extrahieren
data.df <- cbind(data.df, extract(composite_0922.r, xy))
data.df <- na.omit(data.df)
data.df$disturbance <- as.factor(data.df$disturbance)
table(data.df$disturbance) # Zeigt, wie viele Punkte wir pro Klasse haben
summary(data.df)
# Spektrale Indizes als potenzielle Merkmale berechnen
# Normalisierter Differenz-Vegetationsindex (NDVI)
NDVI <- function(x) (x$NIR - x$RED) / (x$NIR + x$RED)
data.df$NDVI <- NDVI(data.df)
# Betrachte die folgenden Plots und Ergebnisse
# --> Diskutiere, welches Merkmal geeignet wäre, um Brandflächen
# und Waldbestände mit abgestorbenen Fichten zu klassifizieren.
boxplot(NIR ~ disturbance, data=data.df)
boxplot(SWIR2 ~ disturbance, data=data.df)
boxplot(NDVI ~ disturbance, data=data.df)
cor(data.df[,9:15]) # Korrelation zwischen Merkmalen
# 4 Random-Forest-Klassifikator trainieren und validieren ------------------------------
# Im Folgenden verwenden wir die Funktion train() aus dem caret-Paket,
# um einen Random Forest zu trainieren, und trainControl(),
# um das Training mit Cross-Validation durchzuführen.
# Cross-Validation-Kontrolle für das Training einrichten
ctrl <- trainControl(
method = "cv", # Cross-Validation verwenden
number = 10, # Anzahl der Folds (k)
summaryFunction = multiClassSummary # Multi-Klassen-Leistungsmetriken verwenden
)
# Training des Random-Forest-Modells durchführen
model <- train(
disturbance ~ BLUE + GREEN + RED + NIR + SWIR1 + SWIR2 + NDVI,
# Klassifiziert die Störung basierend auf den Bändern und spektralen Indizes
data = data.df, # Verwendeter Datensatz
method = "rf", # Random Forest aus dem randomForest-Paket verwenden
trControl = ctrl # Gewählte Kontrollstrategie verwenden
model
confusionMatrix(model, norm="none") # Konfusionsmatrix ausgeben
# 5 Random-Forest-Modell auf das gesamte Untersuchungsgebiet anwenden --------------------
# NDVI für das gesamte Landsat-Bild berechnen
composite_0922.r$NDVI <- NDVI(composite_0922.r)
plot(composite_0922.r$NDVI) # Was siehst du? Was könnte problematisch sein?
# Das trainierte Random-Forest-Modell auf das gesamte Bild anwenden
disturbance.r <- predict(composite_0922.r, model, type="raw", na.rm=TRUE)
plot(disturbance.r, mar=c(2, 2, 2, 7))
# Was hältst du von diesem Klassifikationsergebnis?
# Was hat gut funktioniert? Was nicht?
# 6 CLC-Karte mit den zwei Störungsklassen aktualisieren ------------------------
# Die CLC-Karte mit neuen Klassen für Brandflächen und abgestorbene Fichten aktualisieren
clc_new.r <- resample(clc.r, composite_0922.r, method="near") # CLC-Karte auf Landsat-Auflösung resamplen
# Aktualisierung für die Brandfläche
clc_new.r[clc_new.r == "Woody needle leaved trees" & disturbance.r == "burnt"] <- 91
clc_new.r[clc_new.r == "Woody Broadleaved deciduous trees" & disturbance.r == "burnt"] <- 91
clc_new.r[clc_new.r == "Permanent herbaceous" & disturbance.r == "burnt"] <- 91
clc_new.r[clc_new.r == "Non and sparsely vegetated" & disturbance.r == "burnt"] <- 91
clc_new.r[clc_new.r == "Low-growing woody" & disturbance.r == "burnt"] <- 91
# Aktualisierung für abgestorbene Fichtenbestände
clc_new.r[clc_new.r == "Woody needle leaved trees" & disturbance.r == "deadSpruce"] <- 92
# Neue Legende und Farbtabelle erstellen
clc_new.lgd <- rbind(clc.lgd,
data.frame(i=91, class="Brandfläche", color="#000000"),
data.frame(i=92, class="Abgestorbene Fichte", color="#FF7777"))
levels(clc_new.r) <- clc_new.lgd
coltab(clc_new.r) <- clc_new.lgd[,-2]
plot(clc_new.r, mar=c(2, 2, 2, 9))
# Was hältst du von dieser 'finalen' Karte?
# Was hat gut funktioniert? Was hat nicht gut funktioniert?
# 1 Notwendige Pakete installieren
library(terra) # zum Arbeiten mit Rasterdaten
library(randomForest) # Random-Forest-Klassifikator
library(caret) # für Trainings-, Validierungs- und Regressionsalgorithmen
library(dplyr) # zum Arbeiten und Manipulieren von Data Frames
# 2 Daten einlesen und vorbereiten ---------------------------------------------
# 2.1 Satellitendaten / Karten einlesen -----------------------------------------
# Landsat-9 vom 5. September 2022, plottet: RGB, lineare Kontraststreckung, Titel
# Das Koordinatenreferenzsystem der Landsat-Daten wird für alle weiteren Datensätze verwendet
Aufgaben:
# Verbessern Sie das Ergebnis der Klassifikation und der finalen Karte.
# Führen Sie dazu die folgenden Analysen (7.1 und 7.2) durch:
------------------
7.1
NBR <- function(x) (x$NIR - x$SWIR2) / (x$NIR + x$SWIR2)
data.df$NBR <- NBR(data.df)
—————————————————————————————————-
a)
boxplot(NBR ~ disturbance, data=data.df)
cor(data.df[,9:16]) # Korrelation zwischen den Merkmalen
b)
model2 <- train(
disturbance ~ BLUE + GREEN + RED + NIR + SWIR1 + SWIR2 + NDVI + NBR,
data = data.df,
method = "rf",
trControl = ctrl
model2
confusionMatrix(model2, norm="none")
—————————————————————————————————
c)
print(model$results)
print(model2$results)
composite_0922.r$NBR <- NBR(composite_0922.r)
plot(composite_0922.r$NBR)
disturbance2.r <- predict(composite_0922.r, model2, type="raw", na.rm=TRUE)
plot(disturbance2.r, mar=c(2, 2, 2, 7))
d)
clc_new2.r <- resample(clc.r, composite_0922.r, method="near") # resample the CLC map to the Landsat image
clc_new2.r[clc_new2.r == "Woody needle leaved trees" & disturbance2.r == "burnt"] <- 91
clc_new2.r[clc_new2.r == "Woody Broadleaved deciduous trees" & disturbance2.r == "burnt"] <- 91
clc_new2.r[clc_new2.r == "Permanent herbaceous" & disturbance2.r == "burnt"] <- 91
clc_new2.r[clc_new2.r == "Non and sparsely vegetated" & disturbance2.r == "burnt"] <- 91
clc_new2.r[clc_new2.r == "Low-growing woody" & disturbance2.r == "burnt"] <- 91
clc_new2.r[clc_new2.r == "Woody needle leaved trees" & disturbance2.r == "deadSpruce"] <- 92
clc_new2.lgd <- rbind(clc.lgd,
data.frame(i=91, class="Burnt area", color="#000000"),
data.frame(i=92, class="Dead spruce", color="#FF7777"))
levels(clc_new2.r) <- clc_new2.lgd
coltab(clc_new2.r) <- clc_new2.lgd[,-2]
plot(clc_new2.r, mar=c(2, 2, 2, 9))
# Berechnen Sie den Normalized Burn Ratio (NBR) (siehe z. B. Beetz et al. 2024) auf
# die gleiche Weise, wie oben der NDVI berechnet wurde:
# NBR <- function(x) ...
# data.df$NBR <- NBR(data.df)
- - - - - - - - - - -
—> NBR berechnet
# Untersuchen Sie die Korrelation des NBR mit anderen Merkmalen und wie gut sich
# der NBR zur Trennung der Störungsklassen eignet.
— - - - - - - - - - -
~ bedeutet „in Abhängigkeit von“.
~
NIR ~ disturbance zeigt die Verteilung des NIR-Bandes für jede Störungsklasse (disturbance).
NIR ~ disturbance
disturbance
erstellt Boxplot (mit Minimum, Median, Maximum, Ausreißer)
# Trainieren Sie das Random-Forest-Modell erneut und verwenden Sie zusätzlich den NBR.
- - - - - - - - - - - - -
disturbance ~ BLUE + GREEN + RED + NIR + SWIR1 + SWIR2 + NDVI + NBR, # klassifiziert die Störung basierend auf den Bändern und spektralen Indizes
data = data.df, # der zu verwendende Datensatz
method = "rf", # Verwendung von Random Forest aus dem randomForest-Paket
trControl = ctrl # Verwendung der gewählten Kontrollparameter
confusionMatrix(model2, norm="none") #prüfen, wie gut das Modell die Klassen richtig vorhersagt
Ziel: sehen, ob die Aufnahme von NBR die Klassifikation verbessert hat.
#Vergleich der Leistungsmetriken des Modelltrainings mit den vorherigen Ergebnissen (d.h. ohne NBR).
print(model$results) #model = ursprüngliches Modell ohne NBR
print(model2$results) #model2 = neues Modell mit NBR
#print(...$results) zeigt die Leistungskennzahlen des Modells an
print(...$results)
# Anwendung des trainierten Modells auf das gesamte Kompositbild
composite_0922.r$NBR <- NBR(composite_0922.r) #berechnet
plot(composite_0922.r$NBR) #berechnet den Normalized Burn Ratio (NBR) für das gesamte Satellitenbild-Raster composite_0922.r
composite_0922.r
#plot(composite_0922.r$NBR) erzeugt eine Bildschirmdarstellung des neu berechneten NBR-Bands als Rastergrafik. Die Standard-Farbskala zeigt typischerweise hohe Werte (gesunde Vegetation) in Grün-/Gelbtönen und niedrige Werte (potenzielle Störungen/Brand) in Rot-/Dunkeltönen, was eine schnelle visuelle Einschätzung der Störungsverteilung ermöglicht.
#Vorhersage der Störungen:
disturbance2.r <- predict(composite_0922.r, model2, type="raw", na.rm=TRUE) #wendet das trainierte Random-Forest-Modell model2 (das nun NBR enthält) auf das gesamte Raster an.
plot(disturbance2.r, mar=c(2, 2, 2, 7))#rendert die vorhergesagte Störungskarte mit angepassten Randeinstellungen (mardefiniert Abstände: unten 2, links 2, oben 2, rechts 7 Zeilen für eine Legende)
mar
d) # Aktualisiere die CLC-Karte und diskutiere die Ergebnisse.
clc_new2.r <- resample(clc.r, composite_0922.r, method="near")
Funktion: Passt die Auflösung der ursprünglichen CLC-Karte (clc.r) an die des Landsat-Bildes (composite_0922.r) an. method="near" verwendet Nearest Neighbor-Interpolation (schnell, erhält Klassenwerte).
clc.r
method="near"
2. Klassenaktualisierung für Brandflächen (Klasse 91)
clc_new2.r[clc_new2.r == "Woody needle leaved trees" & disturbance2.r == "burnt"]
<- 91 # ... weitere Vegetationsklassen
<- 91
# ... weitere Vegetationsklassen
Funktion: Ersetzt 5 verschiedene Vegetationsklassen (Nadelwald, Laubwald, Kräuter, etc.) durch neue Klasse 91 ("Burnt area"), wo die Störungsklassifikation "burnt" anzeigt.
3. Klassenaktualisierung für abgestorbene Fichten (Klasse 92)
Funktion: Ersetzt Nadelwaldklassen speziell bei "deadSpruce"-Störung durch neue Klasse 92.
4. Neue Legendeneinträge hinzufügen
clc_new2.lgd <- rbind(clc.lgd, data.frame(i=91, class="Burnt area", color="#000000"), data.frame(i=92, class="Dead spruce", color="#FF7777"))
Funktion: Erweitert die bestehende Legende um zwei neue Klassen:
91 = "Burnt area" (Schwarz #000000)
#000000
92 = "Dead spruce" (Rot #FF7777)
#FF7777
5. Legende und Farbtabelle zuweisen
levels(clc_new2.r) <- clc_new2.lgd coltab(clc_new2.r) <- clc_new2.lgd[,-2]
Funktion:
levels() weist die vollständige Legende (inkl. Namen) zu
levels()
coltab() weist Farben zu (ohne Klassennamen-Spalte)
coltab()
6. Finale Visualisierung
Funktion: Zeigt die aktualisierte CLC-Karte mit neuen Brand- (schwarz) und Totholz-Klassen (rot). mar=c(2,2,2,9) gibt extra Platz rechts für die erweiterte Legende.
mar=c(2,2,2,9)
Aufgaben 7.2
Prüfen, ob das Difference Normalised Burn Ratio die Ergebnisse verbessern kann.
Zusätzlich die Landsat‑8- und Sentinel‑2-Bilder vom 10. September 2021 einlesen
und ein Komposit für 2021 erstellen.
#
l8_100921.r <- rast("data/l8_100921.tif")
plotRGB(l8_100921.r, stretch="lin", main="Landsat 8 (2021-09-10)")
s2a_100921.r <- rast("data/s2a100921_clouds.tif")
plotRGB(s2a_100921.r, stretch="lin", main="Sentinel-2a (2021-09-10)")
composite_0921.r <- l8_100921.r
s2a_100921.r <- resample(s2a_100921.r, composite_0921.r)
composite_0921.r$Blue <- mean(c(l8_100921.r$Blue, s2a_100921.r$Blue), na.rm=TRUE)
composite_0921.r$Green <- mean(c(l8_100921.r$Green, s2a_100921.r$Green), na.rm=TRUE)
composite_0921.r$Red <- mean(c(l8_100921.r$Red, s2a_100921.r$Red), na.rm=TRUE)
composite_0921.r$NIR <- mean(c(l8_100921.r$NIR, s2a_100921.r$NIR), na.rm=TRUE)
composite_0921.r$SWIR1 <- mean(c(l8_100921.r$SWIR1, s2a_100921.r$SWIR1), na.rm=TRUE)
composite_0921.r$SWIR2 <- mean(c(l8_100921.r$SWIR2, s2a_100921.r$SWIR2), na.rm=TRUE)
plot(composite_0921.r)
plotRGB(composite_0921.r, stretch="lin", main="L8 + S2 composite (2021-09)")
————————————————————————————————-
composite_0921.r$NBR <- NBR(composite_0921.r)
plot(composite_0921.r$NBR)
colnames(data.df)[which(colnames(data.df) == 'NBR')] <- 'NBRpost'
data.df$NBRpre <- extract(composite_0921.r$NBR, xy[data.df$X])[,1]
colnames(data.df)[which(colnames(data.df) == 'NBR')] <- 'NBRpre'
data.df$dNBR <- data.df$NBRpre - data.df$NBRpost
data.df$NBRpre <- NULL
data.df$NBRpost <- NULL
boxplot(dNBR ~ disturbance, data=data.df)
cor(data.df[,9:16]) # correlation between features
model3 <- train(
disturbance ~ BLUE + GREEN + RED + NIR + SWIR1 + SWIR2 + NDVI + dNBR,
# classifies the disturbance based on the bands and spectral indices
data = data.df, # our dataset to be used
method = "rf", # use random forest from the randomForest package
trControl = ctrl # use the chosen control setup
model3
confusionMatrix(model3, norm="none")
print(model3$results)
composite_0922.r$dNBR <- composite_0921.r$NBR - composite_0922.r$NBR
plot(composite_0922.r$dNBR, range=c(-1,1))
disturbance3.r <- predict(composite_0922.r, model3, type="raw", na.rm=TRUE)
plot(disturbance3.r, mar=c(2, 2, 2, 7))
———————————————————————————————-
clc_new3.r <- resample(clc.r, composite_0922.r, method="near") # resample the CLC map to the Landsat image
clc_new3.r[clc_new3.r == "Woody needle leaved trees" & disturbance3.r == "burnt"] <- 91
clc_new3.r[clc_new3.r == "Woody Broadleaved deciduous trees" & disturbance3.r == "burnt"] <- 91
clc_new3.r[clc_new3.r == "Permanent herbaceous" & disturbance3.r == "burnt"] <- 91
clc_new3.r[clc_new3.r == "Non and sparsely vegetated" & disturbance3.r == "burnt"] <- 91
clc_new3.r[clc_new3.r == "Low-growing woody" & disturbance3.r == "burnt"] <- 91
clc_new3.r[clc_new3.r == "Woody needle leaved trees" & disturbance3.r == "deadSpruce"] <- 92
clc_new3.lgd <- rbind(clc.lgd,
levels(clc_new3.r) <- clc_new3.lgd
coltab(clc_new3.r) <- clc_new3.lgd[,-2]
plot(clc_new3.r, mar=c(2, 2, 2, 9))
Berechnet den NBR aus einem Multiband-Raster.
Extrahiert Rasterwerte für Punktdaten.
Berechnet dNBR zur Quantifizierung von Störungen.
Bereinigt den Datensatz.
Vergleicht Spektralbänder mittels Boxplots.
Analysiert Korrelationen zwischen spektralen Merkmalen.
# Landsat-8 vom 10. September 2021
#rast() lädt eine Rasterdatei (GeoTIFF) in R. Hier wird ein Landsat 8 Bild vom 10. September 2021 geladen.
rast()
plotRGB() zeigt das Bild als Farbkomposit (RGB), stretch="lin" sorgt für lineare Kontrasterhöhung für bessere Sichtbarkeit.
plotRGB()
stretch="lin"
#Sentinel-2a vom 10. September 2021
#Erstellen eines kombinierten Bildes (Composite)
# Verwende das Landsat-Bild als Grundlage für das Kombinationsbild (Composite)
#composite_0921.r wird als Basis das Landsat-Bild nehmen.
#composite_0921.r
# Passe die Sentinel-2-Daten an die Auflösung und das Raster des Landsat-Bildes an
#resample() passt die Auflösung des Sentinel-2 Bildes an die des Landsat-Bildes an, damit die Pixel direkt verglichen werden können.
resample()
#Layer zusammenführen: Ergebnis ein Landsat + Sentinel-2 Composite vom September 2021.
#Die Bänder Blue, Green, Red, NIR, SWIR1, SWIR2 werden pixelweise gemittelt, um ein kombiniertes Bild zu erzeugen.
Blue
Green
Red
NIR
SWIR1
SWIR2
#na.rm=TRUE ignoriert fehlende Werte.
#na.rm=TRUE
Berechne das differenzierte normalisierte Brandverhältnis (dNBR) aus den NBR-Werten von 2021 und 2022 (siehe z. B. Beetz et al. 2024 für Details).
plot(composite_0921.r$NBR) #NBR() berechnet den Normalized Burn Ratio (NBR).
NBR()
composite_0921.r ist vermutlich ein Multiband-Raster.
composite_0921.r
Mit $NBR wird eine neue Raster-Band-Ebene namens "NBR" hinzugefügt. Das Ergebnis wird im Rasterobjekt gespeichert.
$NBR
"NBR"
plot() visualisiert die neue NBR-Rasterebene.
plot()
👉 „Die Spalte NBR im Dataframe wird in NBRpost umbenannt.“
colnames(data.df) gibt alle Spaltennamen zurück.
colnames(data.df)
which(...) sucht die Position der Spalte "NBR".
which(...)
Diese wird in "NBRpost" umbenannt.
"NBRpost"
👉 „Für die angegebenen Koordinaten werden NBR-Werte aus dem Raster extrahiert und als NBRpre im Dataframe gespeichert.“
extract() liest Rasterwerte an bestimmten Koordinaten aus.
extract()
xy[data.df$X] enthält vermutlich Punktkoordinaten.
xy[data.df$X]
[,1] wählt die erste Spalte des Ergebnisses.
[,1]
Die Werte werden als neue Spalte NBRpre gespeichert.
NBRpre
colnames(data.df)[which(colnames(data.df) == 'NBR')] <- 'NBRpre
Hier wird nochmal eine Spalte "NBR" in "NBRpre" umbenannt.
"NBRpre"
Entfernt alle Zeilen mit NA-Werten.
NA
👉 „Alle Beobachtungen mit fehlenden Werten werden aus dem Datensatz entfernt.“
Differenz zwischen Vorher- und Nachher-NBR.
dNBR wird oft zur Brandintensitätsanalyse verwendet.
👉 „Die Differenz zwischen NBRpre und NBRpost wird berechnet und als dNBR gespeichert.“
Entfernt beide Spalten aus dem Dataframe.
👉 „Die Zwischenvariablen NBRpre und NBRpost werden aus dem Datensatz gelöscht.“
Formel-Schreibweise: Variable ~ Gruppierungsvariable
Zeigt die Verteilung von NIR getrennt nach disturbance-Klassen.
👉 „Es wird ein Boxplot erstellt, der die Verteilung von NIR-Werten nach Störungskategorie zeigt.“
→ Vergleich der Brandintensität zwischen Störungsklassen.
cor(data.df[,9:16])
[,9:16] wählt Spalten 9 bis 16.
cor() berechnet die Pearson-Korrelation.
Ergebnis ist eine Korrelationsmatrix.
👉 „Es wird die lineare Korrelation zwischen den Merkmalen in den Spalten 9 bis 16 berechnet.“
Wiederhole anschließend die Analyse aus Abschnitt 7.1, betrachte die Modelltrainings‑Performance und untersuche die finale Karte. Hat die Verwendung des dNBR die Ergebnisse verbessert?
train()
Die Funktion stammt aus dem Paket caret. Sie trainiert ein Machine-Learning-Modell.
caret
disturbance ~ BLUE + GREEN + RED + NIR + SWIR1 + SWIR2 + NDVI + dNBR
👉 Bedeutung:
disturbance = Zielvariable (abhängige Variable)
Alle anderen = Prädiktoren (Features)
In der Klausur:
"Es wird ein Klassifikationsmodell trainiert, das die Störungsklasse anhand spektraler Bänder und Indizes vorhersagt.“
data = data.df
Der Datensatz, der für das Training verwendet wird.
method = "rf"
"rf" steht für Random Forest
"rf"
Ensemble-Methode aus vielen Entscheidungsbäumen
Reduziert Overfitting durch Bagging
„Es wird ein Random-Forest-Klassifikator verwendet.“
ctrl enthält die Trainingsstrategie
ctrl
z. B. Cross-Validation oder Bootstrapping
steuert Resampling und Modellbewertung
„Die Modellbewertung erfolgt mittels der definierten Cross-Validation-Einstellungen.“
👉 Gibt Trainingsmetriken aus, z. B.:
Accuracy
Kappa
Resampling-Ergebnisse
„Die Modellleistung auf Basis der Resampling-Ergebnisse wird angezeigt.“
Die Konfusionsmatrix zeigt:
True Positives
False Positives
True Negatives
False Negatives
norm="none" bedeutet:
→ Keine Normalisierung, absolute Zahlen
In der Klausur:Die Konfusionsmatrix zeigt die Klassifikationsleistung des Modells anhand der korrekt und falsch klassifizierten Beobachtungen.“
Vergleiche die Leistungskennzahlen des Modelltrainings mit den vorherigen Ergebnissen (also ohne dNBR) und beschreibe, ob die Verwendung von dNBR den Klassifikator verbessert hat.
Es werden die Trainings- und Validierungsergebnisse der drei Modelle ausgegeben, um ihre Performance zu vergleichen.“
Wende das trainierte Modell auf das vollständige Kompositbild an.
🔎 Erklärung
Es wird die Differenz zwischen zwei NBR-Rastern berechnet.
composite_0921.r$NBR = NBR vor Störung
composite_0922.r$NBR = NBR nach Störung
Ergebnis = dNBR-Raster
Formel:
dNBR=NBRpre−NBRpost
dNBR=NBRpre−NBRpost
👉 Interpretation:
Hohe positive Werte → starke Vegetationsverluste
Werte nahe 0 → keine Veränderung
Negative Werte → Vegetationszunahme
📌 Klausurformulierung:
„Es wird ein Differenz-NBR-Raster erzeugt, das spektrale Veränderungen zwischen zwei Zeitpunkten quantifiziert.“
Aktualisiere die CLC-Karte und diskutiere die Ergebnisse.
clc_new3.r <- resample(clc.r, composite_0922.r, method="near")
clc.r = ursprüngliche CORINE Land Cover (CLC)-Karte
composite_0922.r = Landsat-Raster
resample() passt:
Auflösung
Rastergröße
Pixelraster
an das Zielraster an.
🔹 method="near"
Nearest-Neighbour-Interpolation
Wichtig für kategoriale Daten
Verhindert Mischwerte
„Die CLC-Karte wird mittels Nearest-Neighbour-Resampling an die räumliche Auflösung und Geometrie des Landsat-Rasters angepasst.“
Hier passiert logische Rastermaskierung:
Wenn:
Landbedeckung = bestimmter Vegetationstyp
UND Klassifikation = „burnt“
Dann:
Neue Klasse = 91
👉 Das bedeutet:
Bestimmte Landbedeckung wird überschrieben, wenn sie als verbrannt klassifiziert wurde.
Dasselbe wird für mehrere Vegetationstypen wiederholt.
📌 Klausur:
„Vegetationsklassen werden überschrieben, wenn das Random-Forest-Modell eine Störung (z. B. Brand) erkannt hat.“
👉 Spezifische Baumart wird als neue Klasse definiert.
🔥 Wichtiges Konzept hier
Das ist eine kombinierte Klassifikation:
Basis = Landbedeckungskarte
Zusatzinformation = Störungskarte
Ergebnis = aktualisierte Landnutzungskarte
rbind() fügt neue Zeilen zur bestehenden Legende hinzu Neue Klassen:
91 = Burnt area
92 = Dead spruce
Mit Farbdefinition
„Die bestehende Legendentabelle wird um neue Störungsklassen erweitert.“
Weist Rasterklassen ihre Bezeichnungen zu.
Weist Farbwerte zu.
[,-2] entfernt Spalte mit Klassennamen.
„Den Rasterwerten werden Klassenbezeichnungen und Farbwerte zugewiesen.“
Visualisierung der finalen Landbedeckungskarte.
plot(clc_new.r, ...)
plot(clc_new2.r, ...)
plot(clc_new3.r, ...)
👉 Ziel:
Vergleich verschiedener Modellvarianten:
ohne Indizes
mit NDVI
mit NDVI + dNBR
Übung 5 (folgt)
1. Wählen Sie fünf Modellgütemaße aus Tabelle 1 in Janssen und Heuberger (1995), um verschiedene Eigenschaften der Modellleistung zu quantifizieren.
2. Berechnen Sie auf Grundlage der AGBD-Karte (agbd_pre in t/ha) die gesamte AGBD (in Megatonnen, Mt) der gesamten Untersuchungsregion, die durch die LiDAR-Daten abgedeckt ist. Gehen Sie dabei wie folgt vor: a) Wandeln Sie t/ha in t um. b) Verwenden Sie die Funktion global(), um die gesamte AGBD für die gesamte Region zu berechnen. c) Wandeln Sie t in Mt um.
global()
Wie viel Kohlenstoff ist in der oberirdischen Biomasse gespeichert?
3. Verwenden Sie fun="max" in Schritt 4, um die LiDAR-Höhenkarte zu aggregieren, und vergleichen Sie die geschätzte AGBD.
fun="max"
Last changed9 days ago
