84 lines
2.8 KiB
R
84 lines
2.8 KiB
R
#### 1. Vorbereitung ####
|
||
install.packages('RStoolbox')
|
||
install.packages('raster')
|
||
library(raster)
|
||
library(RStoolbox)
|
||
setwd("~/OneDrive/Dokumente/Fernerkundung/PjS/06_Thermal/Daten Tag")
|
||
|
||
mtl <- 'LC08_L1TP_193023_20190726_20190801_01_T1_MTL.txt'
|
||
metaData <- readMeta(mtl)
|
||
|
||
ls8_2020 <- stack(metaData$DATA$FILES[c(1:7)])
|
||
names(ls8_2020) <- c("ultra_blue", "blue", "green", "red", "NIR", "SWIR1", "SWIR2")
|
||
|
||
#metaData$DATA$FILES
|
||
|
||
ls8_2020_thermal <- stack(metaData$DATA$FILES[9])
|
||
names(ls8_2020_thermal) <- c('TIR1')
|
||
|
||
setwd("~/OneDrive/Dokumente/Fernerkundung/PjS/06_Thermal/Berlin Shapefile")
|
||
BB_shape <- shapefile('Berlin_Bezirke.shp')
|
||
|
||
#### 2. Vorverarbeitung ####
|
||
##### 2.1 Top-Of-Atmosphere Reflectance #####
|
||
offset <- metaData$CALREF$offset[1:7]
|
||
gain <- metaData$CALREF$gain[1:7]
|
||
ls8_2020_toa1 <- gain * ls8_2020 + offset
|
||
|
||
sun_elev_day <- metaData$SOLAR_PARAMETERS['elevation']
|
||
sun_earth_distance_day <- metaData$SOLAR_PARAMETERS['distance']
|
||
ls8_2020_toa <- ls8_2020_toa1 * sun_earth_distance_day^2 / sin(sun_elev_day*pi/180)
|
||
#ls8_2020_toa
|
||
|
||
|
||
|
||
##### 2.2 Brightness Temperature #####
|
||
offset <- metaData$CALRAD$offset[10]
|
||
gain <- metaData$CALRAD$gain[10]
|
||
ls8_2020_thermal.RAD <- gain * ls8_2020_thermal + offset
|
||
|
||
K1 <- metaData$CALBT$K1[1]
|
||
K2 <- metaData$CALBT$K2[1]
|
||
ls8_2020_thermal.BT <- K2/log(K1/ls8_2020_thermal.RAD+1)
|
||
#ls8_2020_thermal.BT
|
||
|
||
#### 3. Shapefile ####
|
||
#plotRGB(ls8_2020_toa, r = 4, g = 3, b = 2, stretch = 'hist')
|
||
#plot(BB_shape, lwd = 2, border = 'red', add = TRUE)
|
||
|
||
ls8_2020_toa.BB <- crop(ls8_2020_toa, BB_shape)
|
||
ls8_2020_thermal.BT.BB <- crop(ls8_2020_thermal.BT, BB_shape)
|
||
|
||
#plotRGB(ls8_2020_toa.BB, r = 4, g = 3, b = 2, stretch = 'hist')
|
||
#plot(BB_shape, lwd = 2, border = 'red', add = TRUE)
|
||
|
||
#plot(ls8_2020_thermal.BT.BB, zlim=c(290,320))
|
||
#plot(BB_shape, lwd = 2, border = 'red', add = TRUE)
|
||
|
||
#### 4. Gleichung zur Berechnung der Oberflaechentemperatur ####
|
||
#LST(°C)=Bt/[1+(w∗Bt/p)∗ln(e)]−273.15
|
||
#e=0,017∗Pv+0,963
|
||
#PV=[(NDVI−NDVImin)/(NDVImax−NDVImin)]
|
||
#### 5. Berechnen des NDVI ####
|
||
ls8.ndvi <- spectralIndices(ls8_2020_toa.BB, red = 'red', nir = 'NIR', indices = 'NDVI')
|
||
#plot(ls8.ndvi)
|
||
#ls8.ndvi
|
||
|
||
#### 6. Berechnen des PV ####
|
||
NDVImin <- minValue(ls8.ndvi)
|
||
NDVImax <- maxValue(ls8.ndvi)
|
||
ls8.PV <- ((ls8.ndvi-NDVImin)/(NDVImax-NDVImin))^2
|
||
#plot(ls8.PV)
|
||
|
||
#### 7. Berechnung des Emissionsgrades ####
|
||
ls8.emissiv <- 0.017 * ls8.PV + 0.963
|
||
plot(ls8.emissiv)
|
||
#### 8. Berechnung der Oberflächentemperatur (LST) ####
|
||
ls8.LST <- (ls8_2020_thermal.BT.BB / (1+(10.8 * ls8_2020_thermal.BT.BB/14388) * log(ls8.emissiv))) - 273.15
|
||
#ls8.LST
|
||
#plot(ls8.LST, zlim=c(20,40))
|
||
#plot(BB_shape, lwd = 2, border = "red", add = TRUE)
|
||
#### 9. Speichern des Ergebnisses ####
|
||
setwd("~/OneDrive/Dokumente/Fernerkundung/PjS/06_Thermal")
|
||
writeRaster(ls8.LST$layer, "ls8_LST.sdat", format = "SAGA", overwrite = TRUE)
|