Test
Tom Ihro
2023-11-14
library(terra)## terra 1.7.55img = rast("Data/sen2_20200821_sen2b.tif")
img = img/100 #setting for reflectivity
setMinMax(img)
#set negative values to nan
#img[img < 0] = NaN
plotRGB(img,4,3,2,stretch = "lin")
landcover <- vect("Data/landcoverNew.kml")
landcover = project(landcover, img)
writeVector(landcover,"Data/landcoverclasses.gpkg",overwrite = TRUE)
plotRGB(img, r=8, g=3, b=2, stretch="lin")
plot(landcover, add=TRUE, col = "green")
img.extract <- extract(img, landcover, fun=mean)
names(img.extract) <- c("ID", "blue", "green", "red", "RE1", "RE2", "RE3", "NIR_b", "NIR_a", "SWIR1", "SWIR2")
row.names(img.extract) <- c(landcover$Name)
img.extract## ID blue green red RE1 RE2
## forest 1 2.0582358 3.612598 2.4310372 5.5804610 17.8737411
## Lake 2 1.9201988 3.096083 1.6392843 2.3779722 0.7482107
## lake2 3 0.9049153 1.151864 0.4000000 0.3567797 0.1218644
## industry 4 16.7102030 18.770074 20.3110701 22.8350738 25.1807011
## Grassland 5 5.8490663 8.552923 11.6102030 14.2185521 17.8430853
## field_barren 6 7.5787073 11.314041 16.0715854 18.6497317 22.7484553
## river_spree 7 0.9058659 1.308436 0.9789385 1.2512849 0.7974302
## residential 8 3.7400000 5.494245 5.3334906 8.8312264 18.7901887
## forrest_tiergarten 9 1.8290234 3.406508 2.3429748 5.6061026 21.4920657
## field_vegetation 10 2.8390579 4.923761 3.4983531 7.0870326 26.2318249
## RE3 NIR_b NIR_a SWIR1 SWIR2
## forest 22.3384397 24.3681738 25.1645301 12.3186259 5.0793528
## Lake 0.7198608 0.6139364 0.3645129 0.2724453 0.2596024
## lake2 0.2227119 0.3532203 0.2138983 0.4489831 0.3301695
## industry 25.1085055 24.2507196 25.9204059 22.8854982 16.0696310
## Grassland 19.6450609 21.7923139 23.1670907 29.6725034 18.0408525
## field_barren 24.8395285 27.1888049 29.0604309 34.3070813 22.3930732
## river_spree 0.7825140 0.6149721 0.8234078 1.8129050 1.6875978
## residential 22.4724528 24.0533019 25.1766981 16.3483962 10.2134906
## forrest_tiergarten 28.0224212 30.0106976 31.2995995 14.6867732 6.4592304
## field_vegetation 33.3780415 34.7781306 35.7589688 16.6397923 8.5873442img.extract[1,2:11]## blue green red RE1 RE2 RE3 NIR_b NIR_a
## forest 2.058236 3.612598 2.431037 5.580461 17.87374 22.33844 24.36817 25.16453
## SWIR1 SWIR2
## forest 12.31863 5.079353lambda<-c(490,560, 665, 705, 740, 783, 842, 865, 1610, 2190)
landcover_col <- c("forest green",
"blue1",
"light blue",
"grey",
"light green",
"beige",
"steelblue",
"red",
"darkseagreen",
"green3")
plot(lambda,img.extract[1,2:11],
col = "forest green",
xlab="Wavelength [nm]",
ylab="reflectance (%)",
main="reflectance spectra",
lwd=2, type='b',
xlim=c(400,2700),
ylim=c(0,40))
#loop through all landcover classes
for(i in 2:10){
lines(lambda, img.extract[i,2:11], col = landcover_col[i],lwd=2, type='b')
}
legend('topright', title='Landcover',
legend=row.names(img.extract),
fill=landcover_col,
cex = 0.8, bty="n")
plot(img[[2]], img[[3]], xlab="green", ylab="red", xlim=c(0,40), ylim=c(0,40), maxcell = 100000, col= "gray35")
ran = 1:5
for(i in ran){
points(img.extract[i,"green"],img.extract[i,"red"], col=landcover_col[i], cex=2, pch=19)
}
legend('topright', title='Landcover',
legend=row.names(img.extract)[ran],
fill=landcover_col[ran],
cex = 0.8, bty="n") 
Homework
1. Feature Space
plot(img[[3]], img[[8]], xlab="red", ylab="NIR", xlim=c(0,40), ylim=c(0,40), maxcell = 100000, col= "gray35")
ran = 1:5 #Define the range of which classes we are looking at
for(i in ran){
points(img.extract[i,"red"],img.extract[i,"NIR_a"], col=landcover_col[i], cex=2, pch=19)
}
legend('topright', title='Landcover',
legend=row.names(img.extract)[ran],
fill=landcover_col[ran],
cex = 0.8, bty="n") 
plot(img[[8]], img[[9]], xlab="NIR_a", ylab="SWIR1", xlim=c(0,40), ylim=c(0,40), maxcell = 100000, col= "gray35")
ran = 1:5 #Define the range of which classes we are looking at
for(i in ran){
points(img.extract[i,"NIR_a"],img.extract[i,"SWIR1"], col=landcover_col[i], cex=2, pch=19)
}
legend('topright', title='Landcover',
legend=row.names(img.extract)[ran],
fill=landcover_col[ran],
cex = 0.8, bty="n") 
Unlike for the Visual bands for the NIR and SWIR bands the bands seem to not be correlated that much. As one can observe by the differing Feature spaces the higher
NDVI
NDVI = (img[[8]] - img[[3]])/(img[[8]] + img[[3]])
#adopt range -1,1 as that is the range that the NDVI function works on.
plot(NDVI,stretch = "lin",range = c(-1,1))
NDVI.class = NDVI
#Classify to 1 -> Water 2 -> Land
NDVI.class[NDVI < 0] = 1
NDVI.class[NDVI >= 0] = 2
color = data.frame(value = c(1,2), color = c("lightblue","brown"))
plot(NDVI.class,col = color$color,type = "classes",levels = c("Water","Land"))
EVI
G = 2.5
C1 = 6
C2 = 7.5
L = 1
rhoNIR = img$sen2_20200821_sen2b_7/100
rhoRED = img$sen2_20200821_sen2b_3/100
rhoBLUE = img$sen2_20200821_sen2b_1/100
EVI = G*((rhoNIR - rhoRED)/(rhoNIR + C1*rhoRED - C2*rhoBLUE+L))
plot(EVI,range = c(-0.2,1))
PVI
X = c(5,25)
Y = c(8,23)
par(mfrow = c(1,1), pty = "s")
plot(img[[3]], img[[8]], xlab="red", ylab="NIR", xlim=c(0,40), ylim=c(0,40), maxcell = 100000, col= "gray35")
grid()
points(X,Y, col = "red")
a = (Y[2] - Y[1])/(X[2]-X[1])
b = Y[1] - a*X[1]
PVI = (img[[8]] - a*img[[3]]-b)/sqrt(1+a^2)
plot(PVI,range=c(-10,40))
Compare Indexes
plot(NDVI,EVI,maxcell = 100000,
xlim=c(-0.2,1),
ylim=c(-1,1),
xlab ="NDVI",
ylab ="EVI",
)
plot(NDVI,PVI,maxcell = 1000000,
xlim=c(-1,1),
ylim=c(-10,40),
xlab ="NDVI",
ylab ="PVI",
)
plot(EVI,PVI,maxcell = 100000,
xlim=c(-1,1),
ylim=c(-10,40),
xlab ="EVI",
ylab ="PVI",
)
We can see that especially with the EVI vs PVI scatter plot that these two seem to correlate the most. Compared to the NDVI vs PVI which has alot of degeneracy at \(PVI \approx -4\) which means that there can be alot of NDVI values corresponding to that specific PVI value.