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)

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

Introduction

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: 12,672 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 Pixel Based Sampling performed in Google Earth Engine and exported as csv.

d <- read.csv("trainPoints_pilot_1.csv")
d <- as.data.frame(d)

#head (d)

dd <- subset (d, select = -c(system.index, MSK_CLDPRB, MSK_CLDPRB_1, MSK_SNWPRB, MSK_SNWPRB_1,
                            OBJECTID, QA60, QA60_1, SCL, SCL_1, cb, cb_1,
                            cloudScore, cloudScore_1, probability, probability_1,
                            waterVapor, waterVapor_1, .geo))

d <- na.omit(dd)


# Reorder columns
t<-d %>% select(Class2, VH_Summer, VH_Winter, deltaVH, VV_Summer, VV_Winter, deltaVV,  DPOL_Summer, DPOL_Winter, deltaDPOL,  DPSVI_Summer, DPSVI_Winter, deltaDPSVI, blue, blue_1,  green, green_1, red, red_1, re1, re1_1, re2, re2_1, re3, re3_1, nir, nir_1, nir2, nir2_1, swir1,swir1_1, swir2, swir2_1, NDVI, NDVI_1, NDVI705, NDVI705_1, NDWI, NDWI_1, NDWI2,NDWI2_1, NBR, NBR_1, nARI, nARI_1, REIP, REIP_1, IRECI, IRECI_1,
             SWI, TPI250, TPI500, TPI750, TRI, VBF)

``

Mean Spectral Reflectance

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

Sentinel-1

RSvars <- c("VH_Summer", "VH_Winter",
            "VV_Summer", "VV_Winter",
            "DPOL_Summer","DPOL_Winter",
            "DPSVI_Summer", "DPSVI_Winter"
             )
varsList <-RSvars
titleVar <- "Sentinel-1 SAR variable Means"
num <-length(varsList)

msrSub<-msr %>% select(all_of(varsList))
#create a matrix
msrSub <-as.matrix((msrSub))

#Create a land cover color palette
mycolor <- c("red", "darkgreen", "orange", "yellow", "grey", "green", "purple", "blue")
#Create an empty plot
plot(1, ylim = c(-35, 0), 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(msrSub), cex=0.8, col= mycolor, lty= 1, lwd=3, bty = "n", ncol = 3)

Sentinel -2

Bands

### Spectral Profile function
specProfile1 <-function(varsList, titleVar, num){
  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", "orange", "yellow", "grey", "green", "purple", "blue")
  #Create an empty plot
  plot(1, ylim = c(0, 0.4), 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("topleft", rownames(msr), cex=0.8, col= mycolor, lty= 1, lwd=3, bty = "n")
}
## S2 Bands  -Summer
RSvars <- c("blue", "green", "red", "re1", "re2", "re3", "nir", "nir2", "swir1", "swir2"
)
varsList <-RSvars
titleVar <-"Sentinel-2 Band Means - Summer"

specProfile1(varsList, titleVar)

##  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
titleVar <-"Sentinel-2 Band Means - Autumn"

specProfile1(varsList, titleVar)

Optical Indices

## Spectral Profile function
specProfile2 <-function(varsList, titleVar){
  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", "orange", "yellow", "grey", "green", "purple", "blue")
  #Create an empty plot
  plot(1, ylim = c(-1.0, 1.0), xlim = c(1,num), xlab = "Indexes", 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("bottomleft", rownames(msr), cex=0.8, col= mycolor, lty= 1, lwd=3, bty = "n")
}

## Optical indices - Summer

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

specProfile2(varsList, titleVar)

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

varsList <-RSvars
titleVar <-"Sentinel-2 Optical Indexes - Autumn"

specProfile2(varsList, titleVar)

Density Estimation Plots

featurePlot(x = t[,2:55],
  y= as.factor(t$Class2),
  plot = "density",
  labels = c("Reflectance", "Density Distribution"),
  scales = list(x=list(relation="free"), y = list(relation= "free")),
  layout = c(3,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", "orange", "yellow", "grey", "green", "purple", "blue"))),
)

Box Plots

featurePlot(x = t[,2:55],
  y= as.factor(t$Class2),
  plot = "box",
  scales = list( y = list(relation= "free"), x=list(rot=90)),
  layout = c(3,1),

)

Correlation Matrix

#Sentinel 1 Bands
s1Bandst <-t %>% select("Class2", "VH_Summer", "VH_Winter",
            "VV_Summer", "VV_Winter",
            "DPOL_Summer","DPOL_Winter",
            "DPSVI_Summer", "DPSVI_Winter"
             )
bandCorrelationsS1<-cor(s1Bandst[,2:9])
corrplot.mixed(bandCorrelationsS1, lower="number",tl.pos = 'lt', upper = "color", order="hclust" )

#bandCorrelationsS2<-cor(t[,14:49])

bandCorrelationsS2b<-cor(t[,14:33])
corrplot.mixed(bandCorrelationsS2b, lower="number",tl.pos = 'lt', upper = "color")

bandCorrelationsS2i<-cor(t[,34:49])
corrplot.mixed(bandCorrelationsS2i, lower="number",tl.pos = 'lt', upper = "color" )

bandCorrelationsTopo<-cor(t[,50:55])
corrplot.mixed(bandCorrelationsTopo, lower="number",tl.pos = 'lt', upper = "color",order="hclust")

#write.table(bandCorrelationsTopo, file = "bandCor1_Topo.txt", sep = ",")