Principal Component Analysis script in R

pca_script.r

# Principal Component Analysis script, taken from: https://www.youtube.com/watch?v=uvLLkDQbdJA&ab_channel=GeomaticaAmbiental
# install libraries
install.packages("sp")
install.packages("raster")
install.packages("RStoolbox")
install.packages("ggplot2")
install.packages("reshape2")
install.packages("factoextra")

# load libraries
library("sp")
library("raster")
library("RStoolbox")
library("ggplot2")
library("reshape2")
library("factoextra")

# Set your working directory
setwd("your_directory")

# load your image in .tif format. I used a Sentinel 2 image
s2 <- stack("image.tif")

# Rename bans
# I excluded the Aerosol (B1), Water Vapor (B9) and Cirrus (B10) bands
names(s2) <- c("B2", "B3", "B4", "B5", "B6", "B7", "B8", "B8A", "B11", "B12")

# Rransform to data frame. I had to run it twice to work.
dtf <- as.data.frame(s2)
dtf <- as.data.frame(s2)

# Correlation between bands
cor(na.omit(dtf))

# PCA stats, gives you variance of each component
acp <- prcomp(na.omit(dtf))
summary(acp)

# Graph relation between components
# add this argument: axes=c(CP1, CP2), if you want to choose different components
# example, to graph components 1 and 3: fviz_pca_var(acp, axes=c(1, 3), col.var = "contrib")
fviz_pca_var(acp, col.var = "contrib")
pca <- rasterPCA(s2)

# To visualize the PCA composite agains the false color image
par(mfrow=c(1,2))
plotRGB(pca$map,1,2,3, stretch="lin", main="PCA b123")
plotRGB(s2,9,8,1, stretch="lin", main="Sentinel b11,8A,2")

# To save the PCA composite as a .tif
PCA_123 = stack(pca$map$PC1, pca$map$PC2, pca$map$PC3)
writeRaster(PCA_123, "PCA_123.tif", drivername="Gtiff")