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: 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

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("trainPoints_pilot_2.csv")
d <- as.data.frame(d)

## Summary of the number of points of each class

summary <- d %>% group_by(Class2) %>% summarize(n())
kable(summary) %>% kable_styling(bootstrap_options = "striped", font_size = 14)
Class2 n()
ag 251
conifer 5441
decid 4969
grass 701
imperv 10
mixedwood 1356
shrub 428
water 413 
#head (d)

dd <- subset (d, select = -c(system.index, 
                            OBJECTID, .geo))

d <- na.omit(dd)
#Create a data frame to store the mean spectral reflectance
msr <-aggregate (d, 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

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", "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_mean", "VH_Winter_mean",
            "VV_Summer_mean", "VV_Winter_mean",
            "DPOL_Summer_mean","DPOL_Winter_mean",
            "DPSVI_Summer_mean", "DPSVI_Winter_mean"
             )
varsList <-SARmeans
titleVar <- "Sentinel-1 SAR Cluster Means"
ylimits <- c(-35,0)
specProfile(varsList, titleVar, ylimits)

SARsd <- c("VH_Summer_stdDev", "VH_Winter_stdDev",
            "VV_Summer_stdDev", "VV_Winter_stdDev",
            "DPOL_Summer_stdDev","DPOL_Winter_stdDev",
            "DPSVI_Summer_stdDev", "DPSVI_Winter_stdDev"
             )
varsList <-SARsd
titleVar <- "Sentinel-1 SAR Cluster Standard Deviations"
ylimits <-c(0,3.5)
specProfile(varsList, titleVar, ylimits)

Density Estimation Plots

#function densplot
densPlot <-function(varsList){
  
  t <-d %>% select("Class2", (all_of(varsList)))
  num <-length(varsList)+1
  
   featurePlot(x = t[,2:num],
     y= as.factor(t$Class2),
     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", "orange", "yellow", "purple", "Blue", "grey", "darkgreen", "green "))),
     par.settings = list(superpose.line = list(col = c("red", "darkgreen", "orange", "yellow", "grey", "green", "purple", "blue"))),
)
}
SARvars <- c("VH_Summer_mean", "VH_Summer_stdDev", "VH_Winter_mean","VH_Winter_stdDev", "deltaVH_mean", "deltaVH_stdDev",
            "VV_Summer_mean","VV_Summer_stdDev", "VV_Winter_mean","VV_Winter_stdDev", "deltaVV_mean", "deltaVV_stdDev",
            "DPOL_Summer_mean","DPOL_Summer_stdDev","DPOL_Winter_mean","DPOL_Winter_stdDev","deltaDPOL_mean", "deltaDPOL_stdDev",
            "DPSVI_Summer_mean","DPSVI_Summer_stdDev", "DPSVI_Winter_mean","DPSVI_Winter_stdDev", "deltaDPSVI_mean", "deltaDPSVI_stdDev"
             )

varsList <-SARvars
densPlot (varsList)

Box Plots

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

)
}

SARvars <- c("VH_Summer_mean", "VH_Summer_stdDev", "VH_Winter_mean","VH_Winter_stdDev", "deltaVH_mean", "deltaVH_stdDev",
            "VV_Summer_mean","VV_Summer_stdDev", "VV_Winter_mean","VV_Winter_stdDev", "deltaVV_mean", "deltaVV_stdDev",
            "DPOL_Summer_mean","DPOL_Summer_stdDev","DPOL_Winter_mean","DPOL_Winter_stdDev","deltaDPOL_mean", "deltaDPOL_stdDev",
            "DPSVI_Summer_mean","DPSVI_Summer_stdDev", "DPSVI_Winter_mean","DPSVI_Winter_stdDev", "deltaDPSVI_mean", "deltaDPSVI_stdDev"
             )

varsList <-SARvars
boxy (varsList)

Violin Plots

RSvars <-SARvars

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

Correlation Matrix

#function Correlate
Correlate <-function(varsList){
 t <-d %>% select("Class2", (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_mean", "deltaVH_mean", "deltaDPOL_mean", "deltaDPSVI_mean")
Correlate (varsList)

varsList <- c(SARsd, "deltaVV_stdDev", "deltaVH_stdDev", "deltaDPOL_stdDev", "deltaDPSVI_stdDev")
Correlate (varsList)

varsList <-SARvars
t <-d %>% select("Class2", (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", "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_mean", "green_mean", "red_mean", "re1_mean", "re2_mean", "re3_mean", "nir_mean", "nir2_mean", "swir1_mean", "swir2_mean"
)
varsList <-RSvars
titleVar <-"Sentinel-2 Band Means - Summer"
ylimits <- c(0,0.4)

specProfile(varsList, titleVar, ylimits)

## Summer SD
varsList <-c("blue_stdDev", "green_stdDev", "red_stdDev", "re1_stdDev", "re2_stdDev", "re3_stdDev", "nir_stdDev", "nir2_stdDev", "swir1_stdDev", "swir2_stdDev"
)
titleVar <- "Sentinel-2 Summer Bands Within-Cluster Standard Deviations"
ylimits <-c(0,0.05)
specProfile(varsList, titleVar, ylimits)

##  S2 Bands  -Fall
RSvars <- c("blue_1_mean", "green_1_mean", "red_1_mean", "re1_1_mean", "re2_1_mean", "re3_1_mean", "nir_1_mean", "nir2_1_mean", "swir1_1_mean", "swir2_1_mean"
)
varsList <-RSvars
ylimits <- c(0,0.4)
titleVar <-"Sentinel-2 Band Means - Autumn"

specProfile(varsList, titleVar, ylimits)

## Autumn SD
varsList <-c("blue_1_stdDev", "green_1_stdDev", "red_1_stdDev", "re1_1_stdDev", "re2_1_stdDev", "re3_1_stdDev",
             "nir_1_stdDev", "nir2_1_stdDev", "swir1_1_stdDev", "swir2_1_stdDev"
)
titleVar <- "Sentinel-2 Autumn Bands Within-Cluster Standard Deviations"
ylimits <-c(0,0.05)
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", "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_mean","NDVI705_mean", "EVI_mean", "NDWI_mean", "NDWI2_mean", "NBR_mean", "nARI_mean", "IRECI_mean"
            )
varsList <-RSvars
titleVar <-"Sentinel-2 Optical Indexes - Summer"
ylimits <- c(-1.0,1.0)
specProfile(varsList, titleVar, ylimits)

## Summer SD
varsList <-c("NDVI_stdDev","NDVI705_stdDev", "EVI_stdDev", "NDWI_stdDev", "NDWI2_stdDev", "NBR_stdDev", "nARI_stdDev", "IRECI_stdDev"
            )

titleVar <- "Summer Indices Within-Cluster Standard Deviations"
ylimits <-c(0,0.1)
specProfile(varsList, titleVar, ylimits)

## Optical indices - Autumn
RSvars <- c("NDVI_1_mean","NDVI705_1_mean", "EVI_1_mean", "NDWI_1_mean", "NDWI2_1_mean", "NBR_1_mean", "nARI_1_mean", "IRECI_1_mean"
            
)

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

specProfile(varsList, titleVar, ylimits)

## Autumn SD
varsList <-c("NDVI_1_stdDev","NDVI705_1_stdDev", "EVI_1_stdDev", "NDWI_1_stdDev", "NDWI2_1_stdDev", "NBR_1_stdDev", "nARI_1_stdDev", "IRECI_1_stdDev"
            )

titleVar <- "Autumn Indices Within-Cluster Standard Deviations"
ylimits <-c(0,0.1)
specProfile(varsList, titleVar, ylimits)

Density Estimation Plots

Bands - Summer

#function densplot
densPlot <-function(varsList){
  
  t <-d %>% select("Class2", (all_of(varsList)))
  num <-length(varsList)+1
  
   featurePlot(x = t[,2:num],
     y= as.factor(t$Class2),
     plot = "density",
     labels = c("Reflectance / Standard Deviation", "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", "orange", "yellow", "grey", "green", "purple", "blue"))),
)
}

S2varsSum <- c("blue_mean", "blue_stdDev","green_mean", "green_stdDev","red_mean",  "red_stdDev",
            "re1_mean", "re1_stdDev","re2_mean", "re2_stdDev","re3_mean", "re3_stdDev",
            "nir_mean", "nir_stdDev","nir2_mean", "nir2_stdDev",
            "swir1_mean", "swir1_stdDev","swir2_mean", "swir2_stdDev"
            )
varsList <-S2varsSum
densPlot (varsList)

Bands - Autumn

S2varsFall <- c(
            "blue_1_mean", "blue_1_stdDev","green_1_mean", "green_1_stdDev","red_1_mean",  "red_1_stdDev",
            "re1_1_mean", "re1_1_stdDev","re2_1_mean", "re2_1_stdDev","re3_1_mean", "re3_1_stdDev",
            "nir_1_mean", "nir_1_stdDev","nir2_1_mean", "nir2_1_stdDev",
            "swir1_1_mean", "swir1_1_stdDev","swir2_1_mean", "swir2_1_stdDev"
            )
varsList <-S2varsFall
densPlot (varsList)

Indices - Summer

IndvarsSum <- c("NDVI_mean","NDVI_stdDev", "NDVI705_mean", "NDVI705_stdDev","EVI_mean", "EVI_stdDev",
             "NDWI_mean","NDWI_stdDev", "NDWI2_mean", "NDWI2_stdDev","NBR_mean", "NBR_stdDev","nARI_mean","nARI_stdDev","REIP_mean", "REIP_stdDev","IRECI_mean","IRECI_stdDev"
                         )
varsList <-IndvarsSum
densPlot (varsList)

Indices - Autumn

IndvarsFall <- c("NDVI_1_mean","NDVI_1_stdDev","NDVI705_1_mean","NDVI705_1_stdDev", "EVI_1_mean","EVI_1_stdDev",
             "NDWI_1_mean", "NDWI_1_stdDev","NDWI2_1_mean","NDWI2_1_stdDev", "NBR_1_mean", "NBR_1_stdDev","nARI_1_mean","nARI_1_stdDev", "REIP_1_mean", "REIP_1_stdDev", "IRECI_1_mean","IRECI_1_stdDev"
                         )
varsList <-IndvarsFall
densPlot (varsList)

Box Plots

Bands - Summer

#function boxy
boxy <-function(varsList){
  
  t <-d %>% select("Class2", (all_of(varsList)))
  num <-length(varsList)+1
  
  featurePlot(x = t[,2:num],
  y= as.factor(t$Class2),
  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,"_mean", sep = "")]-d[,paste0( b, "_1_mean", 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 - Seasonal Differences

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


for (b in bands) {
  
   d[,paste0("diff", b, sep = "")]<- d[,paste0(b,"_mean", sep = "")]-d[,paste0( b, "_1_mean", 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 = "Class2", y = var, fill = "Class2")) + 
      geom_violin() +
      scale_fill_manual(name = "", values=c('red', 'darkgreen', 'orange', 'yellow', 'grey', "green", "purple", "blue")) +
      xlab("Class2")+
      labs(title = (var))
  )
}

Bands- Autumn

RSvars <-S2varsFall

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

Bands - Summer-Autumn Difference

RSvars <-Diffvars

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

Indices: Autumn

RSvars <-IndvarsFall

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

Indices: Seasonal Difference

RSvars <-Inddiff

for (i in 1:length(RSvars)){
  var <- RSvars[i]
  print(
    ggplot(d, aes_string(x = "Class2", y = var, fill = "Class2")) + 
      geom_violin() +
      xlab("Class2")+
      scale_fill_manual(name = "", values=c('red', 'darkgreen', '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("Class2", (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_mean", "green_mean", "red_mean", "re1_mean", "re2_mean", "re3_mean", "nir_mean", "nir2_mean", "swir1_mean", "swir2_mean")
Correlate2 (S2meansSumm )

# S2 band means autumn
S2meansFall <- c("blue_1_mean", "green_1_mean", "red_1_mean", "re1_1_mean", "re2_1_mean", "re3_1_mean", "nir_1_mean", "nir2_1_mean", "swir1_1_mean", "swir2_1_mean")
Correlate2 (S2meansFall)

# Difference variables
varsList <- Diffvars
Correlate2 (varsList)

#Indices
IndMeansSum <- c("NDVI_mean","NDVI705_mean", "EVI_mean", "NDWI_mean", "NDWI2_mean", "NBR_mean", "nARI_mean", "IRECI_mean"
            )
Correlate (IndMeansSum)

IndMeanFall <- c("NDVI_1_mean","NDVI705_1_mean", "EVI_1_mean", "NDWI_1_mean", "NDWI2_1_mean", "NBR_1_mean", "nARI_1_mean", "IRECI_1_mean")
Correlate (IndMeanFall)

varsList <- Inddiff
Correlate2 (varsList)

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


# varsList <-c(S2vars, Indvars, Diffvars, Inddiff)
# t <-d %>% select("Class2", (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("Class2", (all_of(varsList)))
  num <-length(varsList)+1
  
   featurePlot(x = t[,2:num],
     y= as.factor(t$Class2),
     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", "orange", "yellow", "grey", "green", "purple", "blue"))),
)
}
htvars <- c("SWI_mean", "SWI_stdDev", "TPI250_mean", "TPI250_stdDev", "TPI500_mean", "TPI500_stdDev","TPI750_mean", "TPI750_stdDev", 
            "TRI_mean", "TRI_stdDev", "VBF_mean", "VBF_stdDev", "CanHt_mean", "CanHt_stdDev")
           

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 = "Class2", y = var, fill = "Class2")) + 
      geom_violin() +
      scale_fill_manual(name = "", values=c('red', 'darkgreen', 'orange', 'yellow', 'grey', "green", "purple", "blue")) +
      xlab("Class2")+
      labs(title = (var))
  )
}

Correlation Matrix

varsList <- htvars
Correlate (varsList)

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