Fernerkundung und Datenanalyse (von mir)

by Stina O.

Übung 4

Kartierung von Landbedeckungsänderungen (Störungen) in der Böhmisch-Sächsischen Schweiz

Verbesserung der Klassifikation und der finalen Karte

Ziel ist es, das Ergebnis der Klassifikation sowie die finale Karte zu verbessern. Führen Sie dazu folgende Analysen durch:

7.1 Prüfen, ob der Normalised Burn Ratio (NBR) die Ergebnisse verbessern kann

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.

7.2 Prüfen, ob der Difference Normalised Burn Ratio (dNBR) die Ergebnisse verbessern kann

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

l9_050922.r <- rast("data/l9_050922.tif")

plotRGB(l9_050922.r, stretch="lin", main="Landsat 9 (2022-09-05)")

plot(l9_050922.r)

# Das Koordinatenreferenzsystem der Landsat-Daten wird für alle anderen Datensätze verwendet

crs <- crs(l9_050922.r)

crs # UTM-Zone 33N

# Sentinel-2a vom 5. September 2022

s2a_050922.r <- rast("data/s2a050922_clouds.tif")

plotRGB(s2a_050922.r, stretch="lin", main="Sentinel-2a (2022-09-05)")

# Sentinel-2b vom 7. September 2022

s2b_070922.r <- rast("data/s2b070922_clouds.tif")

plotRGB(s2b_070922.r, stretch="lin", main="Sentinel-2b (2022-09-07)")

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

plot(clc.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)

summary(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

source("install_packages.R")

# Notwendige Bibliotheken laden

library(terra) # zum Arbeiten mit Rasterdaten

library(RColorBrewer) # für Farbpaletten

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

l9_050922.r <- rast("data/l9_050922.tif")

plotRGB(l9_050922.r, stretch="lin", main="Landsat 9 (2022-09-05)")

plot(l9_050922.r)

# Das Koordinatenreferenzsystem der Landsat-Daten wird für alle weiteren Datensätze verwendet

crs <- crs(l9_050922.r)

crs # UTM-Zone 33N

# Sentinel-2a vom 5. September 2022

s2a_050922.r <- rast("data/s2a050922_clouds.tif")

plotRGB(s2a_050922.r, stretch="lin", main="Sentinel-2a (2022-09-05)")

# Sentinel-2b vom 7. September 2022

s2b_070922.r <- rast("data/s2b070922_clouds.tif")

plotRGB(s2b_070922.r, stretch="lin", main="Sentinel-2b (2022-09-07)")

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)

summary(data.df)

—————————————————————————————————-

boxplot(NIR ~ disturbance, data=data.df)

boxplot(SWIR2 ~ disturbance, data=data.df)

boxplot(NBR ~ disturbance, data=data.df)

cor(data.df[,9:16]) # Korrelation zwischen den Merkmalen

—————————————————————————————————-

model2 <- train(

disturbance ~ BLUE + GREEN + RED + NIR + SWIR1 + SWIR2 + NDVI + NBR,

data = data.df,

method = "rf",

trControl = ctrl

)

model2

confusionMatrix(model2, norm="none")

—————————————————————————————————

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

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

—————————————————————————————————-

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

erstellt Boxplot (mit Minimum, Median, Maximum, Ausreißer)

—————————————————————————————————

# Trainieren Sie das Random-Forest-Modell erneut und verwenden Sie zusätzlich den NBR.

- - - - - - - - - - - - -

model2 <- train(

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

)

model2

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

# 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

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

—————————————————————————————————-

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

2. Klassenaktualisierung für Brandflächen (Klasse 91)

clc_new2.r[clc_new2.r == "Woody needle leaved trees" & disturbance2.r == "burnt"]

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

clc_new2.r[clc_new2.r == "Woody needle leaved trees" & disturbance2.r == "deadSpruce"] <- 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)
92 = "Dead spruce" (Rot #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
coltab() weist Farben zu (ohne Klassennamen-Spalte)

6. Finale Visualisierung

plot(clc_new2.r, mar=c(2, 2, 2, 9))

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.

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)

summary(data.df)

—————————————————————————————————-

boxplot(NIR ~ disturbance, data=data.df)

boxplot(SWIR2 ~ disturbance, data=data.df)

boxplot(NBR ~ disturbance, data=data.df)

cor(data.df[,9:16]) # Korrelation zwischen den Merkmalen

—————————————————————————————————-

model2 <- train(

disturbance ~ BLUE + GREEN + RED + NIR + SWIR1 + SWIR2 + NDVI + NBR,

data = data.df,

method = "rf",

trControl = ctrl

)

model2

confusionMatrix(model2, norm="none")

—————————————————————————————————

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

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

—————————————————————————————————-

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

erstellt Boxplot (mit Minimum, Median, Maximum, Ausreißer)

—————————————————————————————————

# Trainieren Sie das Random-Forest-Modell erneut und verwenden Sie zusätzlich den NBR.

- - - - - - - - - - - - -

model2 <- train(

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

)

model2

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

# 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

#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.

—————————————————————————————————-

d) # Aktualisiere die CLC-Karte und diskutiere die Ergebnisse.

clc_new2.r <- resample(clc.r, composite_0922.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

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)

clc_new2.r[clc_new2.r == "Woody needle leaved trees" & disturbance2.r == "deadSpruce"] <- 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)
92 = "Dead spruce" (Rot #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
coltab() weist Farben zu (ohne Klassennamen-Spalte)

6. Finale Visualisierung

plot(clc_new2.r, mar=c(2, 2, 2, 9))

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.

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 <- na.omit(data.df)

data.df$dNBR <- data.df$NBRpre - data.df$NBRpost

data.df$NBRpre <- NULL

data.df$NBRpost <- NULL

boxplot(NIR ~ disturbance, data=data.df)

boxplot(SWIR2 ~ disturbance, data=data.df)

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(model$results)

print(model2$results)

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,

data.frame(i=91, class="Burnt area", color="#000000"),

data.frame(i=92, class="Dead spruce", color="#FF7777"))

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

plot(clc_new.r, mar=c(2, 2, 2, 9))

plot(clc_new2.r, mar=c(2, 2, 2, 9))

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

l8_100921.r <- rast("data/l8_100921.tif")

plotRGB(l8_100921.r, stretch="lin", main="Landsat 8 (2021-09-10)")

#rast() lädt eine Rasterdatei (GeoTIFF) in R. Hier wird ein Landsat 8 Bild vom 10. September 2021 geladen.

plotRGB() zeigt das Bild als Farbkomposit (RGB), stretch="lin" sorgt für lineare Kontrasterhöhung für bessere Sichtbarkeit.

#Sentinel-2a vom 10. September 2021

s2a_100921.r <- rast("data/s2a100921_clouds.tif")

plotRGB(s2a_100921.r, stretch="lin", main="Sentinel-2a (2021-09-10)")

#Erstellen eines kombinierten Bildes (Composite)

# Verwende das Landsat-Bild als Grundlage für das Kombinationsbild (Composite)

composite_0921.r <- l8_100921.r

#composite_0921.r wird als Basis das Landsat-Bild nehmen.

# Passe die Sentinel-2-Daten an die Auflösung und das Raster des Landsat-Bildes an

s2a_100921.r <- resample(s2a_100921.r, composite_0921.r)

#resample() passt die Auflösung des Sentinel-2 Bildes an die des Landsat-Bildes an, damit die Pixel direkt verglichen werden können.

#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.

#na.rm=TRUE ignoriert fehlende Werte.

—————————————————————————————————-

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

composite_0921.r$NBR <- NBR(composite_0921.r)

plot(composite_0921.r$NBR) #NBR() berechnet den Normalized Burn Ratio (NBR).

composite_0921.r ist vermutlich ein Multiband-Raster.

Mit $NBR wird eine neue Raster-Band-Ebene namens "NBR" hinzugefügt. Das Ergebnis wird im Rasterobjekt gespeichert.

plot() visualisiert die neue NBR-Rasterebene.

colnames(data.df)[which(colnames(data.df) == 'NBR')] <- 'NBRpost'

👉 „Die Spalte NBR im Dataframe wird in NBRpost umbenannt.“

colnames(data.df) gibt alle Spaltennamen zurück.
which(...) sucht die Position der Spalte "NBR".
Diese wird in "NBRpost" umbenannt.

data.df$NBRpre <- extract(composite_0921.r$NBR, xy[data.df$X])[,1]

👉 „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.
xy[data.df$X] enthält vermutlich Punktkoordinaten.
[,1] wählt die erste Spalte des Ergebnisses.
Die Werte werden als neue Spalte NBRpre gespeichert.

colnames(data.df)[which(colnames(data.df) == 'NBR')] <- 'NBRpre

Hier wird nochmal eine Spalte "NBR" in "NBRpre" umbenannt.

data.df <- na.omit(data.df)

Entfernt alle Zeilen mit NA-Werten.

👉 „Alle Beobachtungen mit fehlenden Werten werden aus dem Datensatz entfernt.“

data.df$dNBR <- data.df$NBRpre - data.df$NBRpost

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.“

data.df$NBRpre <- NULL

data.df$NBRpost <- NULL

Entfernt beide Spalten aus dem Dataframe.

👉 „Die Zwischenvariablen NBRpre und NBRpost werden aus dem Datensatz gelöscht.“

boxplot(NIR ~ disturbance, data=data.df)

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.“

boxplot(dNBR ~ disturbance, data=data.df)

→ 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?

model3 <- train(

`train()`

Die Funktion stammt aus dem Paket caret. Sie trainiert ein Machine-Learning-Modell.

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
Ensemble-Methode aus vielen Entscheidungsbäumen
Reduziert Overfitting durch Bagging

In der Klausur:

„Es wird ein Random-Forest-Klassifikator verwendet.“

`trControl = ctrl`

ctrl enthält die Trainingsstrategie
z. B. Cross-Validation oder Bootstrapping
steuert Resampling und Modellbewertung

In der Klausur:

„Die Modellbewertung erfolgt mittels der definierten Cross-Validation-Einstellungen.“

model3

👉 Gibt Trainingsmetriken aus, z. B.:

Accuracy

Kappa

Resampling-Ergebnisse

In der Klausur:

„Die Modellleistung auf Basis der Resampling-Ergebnisse wird angezeigt.“

confusionMatrix(model3, norm="none")

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.

print(model$results)

print(model2$results)

print(model3$results)

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.

composite_0922.r$dNBR <- composite_0921.r$NBR - composite_0922.r$NBR

🔎 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

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

🔎 Erklärung

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

📌 Klausurformulierung:

„Die CLC-Karte wird mittels Nearest-Neighbour-Resampling an die räumliche Auflösung und Geometrie des Landsat-Rasters angepasst.“

clc_new3.r[clc_new3.r == "Woody needle leaved trees" & disturbance3.r == "burnt"] <- 91

🔎 Erklärung

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.“

clc_new3.r[clc_new3.r == "Woody needle leaved trees" & disturbance3.r == "deadSpruce"] <- 92

👉 Spezifische Baumart wird als neue Klasse definiert.

🔥 Wichtiges Konzept hier

Das ist eine kombinierte Klassifikation:

Basis = Landbedeckungskarte

Zusatzinformation = Störungskarte

Ergebnis = aktualisierte Landnutzungskarte

clc_new3.lgd <- rbind(clc.lgd,

data.frame(i=91, class="Burnt area", color="#000000"),

data.frame(i=92, class="Dead spruce", color="#FF7777"))

🔎 Erklärung

rbind() fügt neue Zeilen zur bestehenden Legende hinzu Neue Klassen:

91 = Burnt area

92 = Dead spruce

Mit Farbdefinition

📌 Klausur:

„Die bestehende Legendentabelle wird um neue Störungsklassen erweitert.“

levels(clc_new3.r) <- clc_new3.lgd

Weist Rasterklassen ihre Bezeichnungen zu.

coltab(clc_new3.r) <- clc_new3.lgd[,-2]

Weist Farbwerte zu.

[,-2] entfernt Spalte mit Klassennamen.

📌 Klausur:

„Den Rasterwerten werden Klassenbezeichnungen und Farbwerte zugewiesen.“

plot(clc_new3.r, mar=c(2, 2, 2, 9))

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)

Join Course

Preview

Author

Stina O.

Information

Last changed
a day ago

Report course