knitr::opts_chunk$set(echo = FALSE,fig.height = 6, fig.width = 11, fig.align = "center")

library(dplyr)
library(knitr)
library(plotly)
library (sf)
library(caret)
library(corrplot)
library(ggplot2)
library(kableExtra)

setwd("C:/Users/user/OneDrive/ABMI/R_scripts")

Exploratory Data Analysis

Study Area: Boreal Pilot Area, 3.7 million hectares in northern Alberta north of Athabasca, east of Slave Lake, and south of Fort McMurray

Training Data: 8469 random points from interpreted photo 3 by 7 km photo-plots distributed across study area Spectral Data: - Sentinel-1 Summer 2021 and Winter 2020-2021 mean composites and derived indices - Sentinel-2 Summer 2019-2021 and Autumn 2019-2021 mean composites and derived indices - Compiled in Google Earth Engine

Objected Based Sampling Means and standard deviations of clusters created at 5m scale in Google Earth Engine. The clusters (objects) average 2 hectares in size, with a minimum of 25 sq.m and a maximum of 9.8 ha.

d <- read.csv("G:\\My Drive\\BorealPilotArea_UplandClass\\trainPoints_sampled_Mar30_toDrive.csv")
d <- as.data.frame(d)

## Summary of the number of points of each class

summary <- d %>% group_by(Class7) %>% summarize(n())
kable(summary) %>% kable_styling(bootstrap_options = "striped", font_size = 14)
Class7 n()
burn 407
Conifer 2026
crop 461
Decid 1527
grass 883
imperv 102
Mixedwood 2051
shrub 589
water 423 
#head (d)

dd <- subset (d, select = -c(system.index,  Class3, Class5, Class6,
                            .geo))

d <- na.omit(dd)
#Create a data frame to store the mean spectral reflectance
msr <-aggregate (d, list(d$Class7), mean, na.rm = TRUE)
#remove the first column in order to use the actual land cover names
rownames(msr)<-msr[,1]
msr<-msr[,-1]
#head (d)

Sentinel-1

Spectral Profiles

### Spectral Profile function
specProfile <-function(varsList, titleVar, ylimits){
  msrSub<-msr %>% select(all_of(varsList))
  num <-length(varsList)
  #create a matrix
  msrSub <-as.matrix((msrSub))
  
  #Create a land cover color palette
  mycolor <- c("red", "darkgreen", "brown", "orange", "yellow ", "grey", "green" ,"purple", "blue")
  #Create an empty plot
  plot(1, ylim = ylimits, xlim = c(1,num), xlab = "SAR Deriviatives", ylab = "Backscatter", xaxt= 'n')
  #custom x-axis
  axis(1, at=1: num, lab=colnames(msrSub))
  #add spectral reflectance
  for (i in 1: nrow(msrSub)){
    lines(msrSub[i,], type = "o", lwd = 3, lty = 1, col = mycolor[i])
  }
  title(main = paste(titleVar), font.main = 2)
  legend("bottom", rownames(msr), cex=0.8, col= mycolor, lty= 1, lwd=3, bty = "n", ncol = 3)
}

SARmeans <- c("VH_Summer", "VH_Winter",
            "VV_Summer", "VV_Winter",
            "DPOL_Summer","DPOL_Winter",
            "DPSVI_Summer", "DPSVI_Winter"
             )
varsList <-SARmeans
titleVar <- "Sentinel-1 SAR Cluster Means"
ylimits <- c(-35,0)
specProfile(varsList, titleVar, ylimits)

Density Estimation Plots

#function densplot
densPlot <-function(varsList){
  
  t <-d %>% select("Class7", (all_of(varsList)))
  num <-length(varsList)+1
  
   featurePlot(x = t[,2:num],
     y= as.factor(t$Class7),
     plot = "density",
     labels = c("Backscatter / Standard Deviation", "Density Distribution"),
     scales = list(x=list(relation="free"), y = list(relation= "free")),
     layout = c(2,1),
     auto.key = list(columns = 3),
     
     par.settings = list(superpose.line = list(col = c("red", "darkgreen", "brown", "orange", "yellow ", "grey", "green" ,"purple", "blue"))),
)
}
SARvars <- c("VH_Summer",  "VH_Winter", "deltaVH", 
            "VV_Summer", "VV_Winter", "deltaVV", 
            "DPOL_Summer","DPOL_Winter","deltaDPOL", 
            "DPSVI_Summer", "DPSVI_Winter", "deltaDPSVI"
             )

varsList <-SARvars
densPlot (varsList)

Box Plots

#function boxy
boxy <-function(varsList){
  
  t <-d %>% select("Class7", (all_of(varsList)))
  num <-length(varsList)+1
  
  featurePlot(x = t[,2:num],
  y= as.factor(t$Class7),
  plot = "box",
  scales = list( y = list(relation= "free"), x=list(rot=90)),
  layout = c(4,1),

)
}

SARvars <- c("VH_Summer",  "VH_Winter", "deltaVH", 
            "VV_Summer", "VV_Winter", "deltaVV",
            "DPOL_Summer","DPOL_Winter","deltaDPOL", 
            "DPSVI_Summer", "DPSVI_Winter", "deltaDPSVI" 
             )

varsList <-SARvars
boxy (varsList)

Violin Plots

RSvars <-SARvars

for (i in 1:length(RSvars)){
  var <- RSvars[i]
  print(
    ggplot(d, aes_string(x = "Class7", y = var, fill = "Class7")) + 
      geom_violin() +
      scale_fill_manual(name = "", values=c("red", "darkgreen", "brown", "orange", "yellow ", "grey", "green" ,"purple", "blue")) +
      xlab("Class7")+
      labs(title = (var))
  )
}

Correlation Matrix

#function Correlate
Correlate <-function(varsList){
 t <-d %>% select("Class7", (all_of(varsList)))
 num <-length(varsList)+1
 bandCorrelationsS1<-cor(t[,2:num])
 corrplot.mixed(bandCorrelationsS1, lower="number",tl.pos = 'lt', upper = "color", order="hclust" )
}

varsList <- c(SARmeans,  "deltaVV", "deltaVH", "deltaDPOL", "deltaDPSVI")
Correlate (varsList)

varsList <-SARvars
t <-d %>% select("Class7", (all_of(varsList)))
num <-length(varsList)+1
bandCorrelationsS1<-cor(t[,2:num])
write.table(bandCorrelationsS1, file = "bandCorS1all.txt", sep = ",")

Sentinel -2

Spectral Profiles

Bands

### Spectral Profile function
specProfile <-function(varsList, titleVar, ylimits){
  msrSub<-msr %>% select(all_of(varsList))
  num <-length(varsList)
  #create a matrix
  msrSub <-as.matrix((msrSub))
  
  #Create a land cover color palette
  mycolor <- c("red", "darkgreen", "brown", "orange", "yellow ", "grey", "green" ,"purple", "blue")
  #Create an empty plot
  plot(1, ylim = ylimits, xlim = c(1,num), xlab = "Bands", ylab = "Reflectance", xaxt= 'n')
  #custom x-axis
  axis(1, at=1: num, lab=colnames(msrSub))
  #add spectral reflectance
  for (i in 1: nrow(msrSub)){
    lines(msrSub[i,], type = "o", lwd = 3, lty = 1, col = mycolor[i])
  }
  title(main = paste(titleVar), font.main = 2)
  legend("top", rownames(msr), cex=0.8, col= mycolor, lty= 1, lwd=3, bty = "n", ncol = 4)
}

#
## S2 Bands  -Summer Mean
RSvars <- c("blue", "green", "red", "re1", "re2", "re3", "nir", "nir2", "swir1", "swir2"
)
varsList <-RSvars
titleVar <-"Sentinel-2 Band Means - Summer"
ylimits <- c(0,0.4)

specProfile(varsList, titleVar, ylimits)

##  S2 Bands  -Fall
RSvars <- c("blue_1", "green_1", "red_1", "re1_1", "re2_1", "re3_1", "nir_1", "nir2_1", "swir1_1", "swir2_1"
)
varsList <-RSvars
ylimits <- c(0,0.4)
titleVar <-"Sentinel-2 Band Means - Autumn"

specProfile(varsList, titleVar, ylimits)

Optical Indices

### Spectral Profile function
specProfile <-function(varsList, titleVar, ylimits){
  msrSub<-msr %>% select(all_of(varsList))
  num <-length(varsList)
  #create a matrix
  msrSub <-as.matrix((msrSub))
  
  #Create a land cover color palette
  mycolor <- c("red", "darkgreen", "brown", "orange", "yellow ", "grey", "green" ,"purple", "blue")
  #Create an empty plot
  plot(1, ylim = ylimits, xlim = c(1,num), xlab = "Indices", ylab = "DN", xaxt= 'n')
  #custom x-axis
  axis(1, at=1: num, lab=colnames(msrSub))
  #add spectral reflectance
  for (i in 1: nrow(msrSub)){
    lines(msrSub[i,], type = "o", lwd = 3, lty = 1, col = mycolor[i])
  }
  title(main = paste(titleVar), font.main = 2)
  legend("top", rownames(msr), cex=0.8, col= mycolor, lty= 1, lwd=3, bty = "n", ncol = 4)
}
## Optical indices - Summer

RSvars <- c("NDVI","NDVI705", "EVI", "NDWI", "NDWI2", "NBR", "nARI", "IRECI"
            )
varsList <-RSvars
titleVar <-"Sentinel-2 Optical Indexes - Summer"
ylimits <- c(-1.0,1.0)
specProfile(varsList, titleVar, ylimits)

## Optical indices - Autumn
RSvars <- c("NDVI_1","NDVI705_1", "EVI_1", "NDWI_1", "NDWI2_1", "NBR_1", "nARI_1", "IRECI_1"
            
)

varsList <-RSvars
ylimits <- c(-1.0,1.0)
titleVar <-"Sentinel-2 Optical Indexes - Autumn"

specProfile(varsList, titleVar, ylimits)

Density Estimation Plots

Bands - Summer

#function densplot
densPlot <-function(varsList){
  
  t <-d %>% select("Class7", (all_of(varsList)))
  num <-length(varsList)+1
  
   featurePlot(x = t[,2:num],
     y= as.factor(t$Class7),
     plot = "density",
     labels = c("Reflectance ", "Density Distribution"),
     scales = list(x=list(relation="free"), y = list(relation= "free")),
     layout = c(2,1),
     auto.key = list(columns = 3),
     lwd = 2,
     #par.settings = list(superpose.line = list(col = c("red", "orange", "yellow", "purple", "Blue", "grey", "darkgreen", "green "))),
     par.settings = list(superpose.line = list(col = c("red", "darkgreen", "brown", "orange", "yellow ", "grey", "green" ,"purple", "blue"))),
)
}

S2varsSum <- c("blue","green" ,"red", 
            "re1","re2","re3", 
            "nir","nir2", 
            "swir1", "swir2"
            )
varsList <-S2varsSum
densPlot (varsList)

Bands - Autumn

S2varsFall <- c(
            "blue_1", "green_1","red_1", 
            "re1_1", "re2_1","re3_1", 
            "nir_1","nir2_1", 
            "swir1_1","swir2_1" 
            )
varsList <-S2varsFall
densPlot (varsList)

Indices - Summer

IndvarsSum <- c("NDVI", "NDVI705", "EVI", 
             "NDWI", "NDWI2", "NBR", "nARI","REIP", "IRECI"
                         )
varsList <-IndvarsSum
densPlot (varsList)

Indices - Autumn

IndvarsFall <- c("NDVI_1","NDVI705_1","EVI_1",
             "NDWI_1", "NDWI2_1", "NBR_1", "nARI_1", "REIP_1",  "IRECI_1"
                         )
varsList <-IndvarsFall
densPlot (varsList)

Box Plots

Bands - Summer

#function boxy
boxy <-function(varsList){
  
  t <-d %>% select("Class7", (all_of(varsList)))
  num <-length(varsList)+1
  
  featurePlot(x = t[,2:num],
  y= as.factor(t$Class7),
  plot = "box",
  scales = list( y = list(relation= "free"), x=list(rot=90)),
  layout = c(4,1),

)
}


varsList <-S2varsSum
boxy (varsList)

Bands - Autumn

varsList <-S2varsFall
boxy (varsList)

Bands: Summer-Autumn Difference

bands <- list ("blue", "green", "red", "re1", "re2", "re3", "nir", "nir2","swir1", "swir2")



for (b in bands) {

  d[,paste0("diff", b, sep = "")]<- d[,paste0(b,"", sep = "")]-d[,paste0( b, "_1", sep = "")]
}

Diffvars <- c("diffblue", "diffgreen", "diffred", "diffre1", "diffre2", 
            "diffre3", "diffnir", "diffnir2", "diffswir1", "diffswir2")

varsList <-Diffvars
boxy (varsList)

Indices

varsList <- c(IndvarsSum, IndvarsFall)
boxy (varsList)

Indices - Summer-Autumn Difference

bands <- list ("NDVI", "NDVI705", "EVI", "NDWI", "NDWI2", "NBR", "nARI",
               "REIP", "IRECI")


for (b in bands) {
  
   d[,paste0("diff", b, sep = "")]<- d[,paste0(b,"", sep = "")]-d[,paste0( b, "_1", sep = "")]
}

Inddiff <- c("diffNDVI", "diffNDVI705", "diffEVI", "diffNDWI", "diffNDWI2", "diffNBR", "diffnARI",
            "diffREIP", "diffIRECI"
)

varsList <-Inddiff
boxy (varsList)

Violin Plots

Bands - Summer

RSvars <-S2varsSum

for (i in 1:length(RSvars)){
  var <- RSvars[i]
  print(
    ggplot(d, aes_string(x = "Class7", y = var, fill = "Class7")) + 
      geom_violin() +
      scale_fill_manual(name = "", values=c("red", "darkgreen", "brown", "orange", "yellow ", "grey", "green" ,"purple", "blue")) +
      xlab("Class7")+
      labs(title = (var))
  )
}

Bands- Autumn

RSvars <-S2varsFall

for (i in 1:length(RSvars)){
  var <- RSvars[i]
  print(
    ggplot(d, aes_string(x = "Class7", y = var, fill = "Class7")) + 
      geom_violin() +
      scale_fill_manual(name = "", values=c("red", "darkgreen", "brown", "orange", "yellow ", "grey", "green" ,"purple", "blue")) +
      xlab("Class7")+
      labs(title = (var))
  )
}

Bands - Summer-Autumn Difference

RSvars <-Diffvars

for (i in 1:length(RSvars)){
  var <- RSvars[i]
  print(
    ggplot(d, aes_string(x = "Class7", y = var, fill = "Class7")) + 
      geom_violin() +
      xlab("Class7")+
      scale_fill_manual(name = "", values=c("red", "darkgreen", "brown", "orange", "yellow ", "grey", "green" ,"purple", "blue")) +
      labs(title = paste("Seasonal Difference ", toupper(substring(var,first = 5))))
  )
}

Indices: Summer

RSvars <-IndvarsSum

for (i in 1:length(RSvars)){
  var <- RSvars[i]
  print(
    ggplot(d, aes_string(x = "Class7", y = var, fill = "Class7")) + 
      geom_violin() +
      scale_fill_manual(name = "", values=c("red", "darkgreen", "brown", "orange", "yellow ", "grey", "green" ,"purple", "blue")) +
      xlab("Class7")+
      labs(title = (var))
  )
}

Indices: Autumn

RSvars <-IndvarsFall

for (i in 1:length(RSvars)){
  var <- RSvars[i]
  print(
    ggplot(d, aes_string(x = "Class7", y = var, fill = "Class7")) + 
      geom_violin() +
      scale_fill_manual(name = "", values=c("red", "darkgreen", "brown", "orange", "yellow ", "grey", "green" ,"purple", "blue")) +
      xlab("Class7")+
      labs(title = (var))
  )
}

Indices: Summer-Autumn Difference

RSvars <-Inddiff

for (i in 1:length(RSvars)){
  var <- RSvars[i]
  print(
    ggplot(d, aes_string(x = "Class7", y = var, fill = "Class7")) + 
      geom_violin() +
      xlab("Class7")+
      scale_fill_manual(name = "", values=c("red", "darkgreen", "brown", "orange", "yellow ", "grey", "green" ,"purple", "blue")) +
      labs(title = paste("Seasonal Difference ", substring(var,first = 5)))
  )
}

Correlation Matrices

#function Correlate
Correlate2 <-function(varsList){
 t <-d %>% select("Class7", (all_of(varsList)))
 num <-length(varsList)+1
 bandCorrelationsS2<-cor(t[,2:num])
 corrplot.mixed(bandCorrelationsS2, lower="number",tl.pos = 'lt', upper = "color" )
}

# S2 band means summer
S2meansSumm  <- c("blue", "green", "red", "re1", "re2", "re3", "nir", "nir2", "swir1", "swir2")
Correlate2 (S2meansSumm )

# S2 band means autumn
S2meansFall <- c("blue_1", "green_1", "red_1", "re1_1", "re2_1", "re3_1", "nir_1", "nir2_1", "swir1_1", "swir2_1")
Correlate2 (S2meansFall)

# Difference variables
varsList <- Diffvars
Correlate2 (varsList)

#Indices
IndMeansSum <- c("NDVI","NDVI705", "EVI", "NDWI", "NDWI2", "NBR", "nARI", "IRECI"
            )
Correlate (IndMeansSum)

IndMeanFall <- c("NDVI_1","NDVI705_1", "EVI_1", "NDWI_1", "NDWI2_1", "NBR_1", "nARI_1", "IRECI_1")
Correlate (IndMeanFall)

varsList <- Inddiff
Correlate2 (varsList)

 # varsList  <- c(Indvars, Inddiff)
 # Correlate2 (varsList)


# varsList <-c(S2vars, Indvars, Diffvars, Inddiff)
# t <-d %>% select("Class7", (all_of(varsList)))
# num <-length(varsList)+1
# bandCorrelationsS2all<-cor(t[,2:num])
# write.table(bandCorrelationsS2all, file = "bandCorS2all.txt", sep = ",")

Topo and Canopy Height

Density Estimation Plots

#function densplot
densPlot <-function(varsList){
  
  t <-d %>% select("Class7", (all_of(varsList)))
  num <-length(varsList)+1
  
   featurePlot(x = t[,2:num],
     y= as.factor(t$Class7),
     plot = "density",
     labels = c("", "Density Distribution"),
     scales = list(x=list(relation="free"), y = list(relation= "free")),
     layout = c(2,1),
     auto.key = list(columns = 3),
     #par.settings = list(superpose.line = list(col = c("red", "orange", "yellow", "purple", "Blue", "grey", "darkgreen", "green "))),
     par.settings = list(superpose.line = list(col = c("red", "darkgreen", "brown", "orange", "yellow ", "grey", "green" ,"purple", "blue"))),
)
}
htvars <- c("SWI",  "VBF",  "CanHt")
           

varsList <-htvars
densPlot (varsList)

Box Plots

varsList <-htvars
boxy (varsList)

Violin Plots

RSvars <-htvars

for (i in 1:length(RSvars)){
  var <- RSvars[i]
  print(
    ggplot(d, aes_string(x = "Class7", y = var, fill = "Class7")) + 
      geom_violin() +
      scale_fill_manual(name = "", values=c("red", "darkgreen", "brown", "orange", "yellow ", "grey", "green" ,"purple", "blue")) +
      xlab("Class7")+
      labs(title = (var))
  )
}

Correlation Matrix

varsList <- htvars
Correlate (varsList)

# varsList <-htvars
# t <-d %>% select("Class7", (all_of(varsList)))
# num <-length(varsList)+1
# bandCorrelations<-cor(t[,2:num])
# write.table(bandCorrelations, file = "bandCorHtall.txt", sep = ",")