A tool for estimating snow cover from alpine meadows using high resolution imagery and random forest

Load Libraries and Images

## Libraries
library(raster)
library(rgdal)
library(rgeos)
library(randomForest)
library(dismo)
library(tidyverse)

## load functions
source("functions.r")

## Load snow image
snow <- brick("images//JumpingPound2012.tif")
plotRGB(snow, stretch = "lin")

## Load meadows
meadows <- readOGR("data//shapeFile//meadowOutline.shp")

## OGR data source with driver: ESRI Shapefile 
## Source: "D:\RStudio\AlpineSnowClassifier\data\shapeFile\meadowOutline.shp", layer: "meadowOutline"
## with 20 features
## It has 2 fields
## Integer64 fields read as strings:  id

## reproject meadows to be same CRS as raster
meadows <- spTransform(meadows, crs(snow))


## create butter around meadows to capture entire area
# plot(meadows)
buffMeadow <- gBuffer(meadows, byid=T,  width=10)
plot(buffMeadow)
plot(meadows, add=T)

## Convert meadow names to characters
buffMeadow$meadow <- as.character(buffMeadow$meadow)
buffMeadow$meadow[buffMeadow$meadow=="g"] <- "littleg"
meadowList <- buffMeadow$meadow
meadowList <- meadowList[!meadowList %in% c("E","D","C")] ## drop Lusk

Example of classifying an individual pixel as snow, forest, rock/bare ground

## Extract a meadow
i=1 ## meadow Z
meadowList[i]

## [1] "Z"

snowMeadow <- crop(snow, buffMeadow[buffMeadow$meadow == meadowList[i],])
snowMeadow <- mask(snowMeadow, buffMeadow[buffMeadow$meadow == meadowList[i],])



## Place 100 random samples to classify
set.seed(11)
spTrain <- sampleRandom(snowMeadow, 100, sp=T)

## Plot locations of pixel
plotRGB(snowMeadow, stretch = "lin")
text(x=coordinates(spTrain)[,1], y=coordinates(spTrain)[,2], labels=1:100)

## Zoom in on a single pixel to classify
j = 1
range = 30 ## range of plot in metres

plotRGB(crop(snowMeadow, extent(
  coordinates(spTrain)[j,1]-range,  ## xmin
  coordinates(spTrain)[j,1]+range,  ## xmax
  coordinates(spTrain)[j,2]-range,  ## ymin
  coordinates(spTrain)[j,2]+range)))  ## ymax
plot(spTrain[j,], add=T)

## Create datatable based on classifications

Apply classified data using Random Forest

### Load classified data
meadowFiles <- list.files("data//snowClass//", pattern=".csv", full.names=T)
trainData <- data.frame()
for(i in 1:length(meadowFiles)){
  temp <- read.table(meadowFiles[i],  header=T, sep=",")
  trainData <- rbind(trainData, temp)
}

## Conduct random forest on RGB values
rf1 <- randomForest(class ~ JumpingPound2012.1 + JumpingPound2012.2 + JumpingPound2012.3, data=trainData)
rf1

## 
## Call:
##  randomForest(formula = class ~ JumpingPound2012.1 + JumpingPound2012.2 +      JumpingPound2012.3, data = trainData) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 1
## 
##         OOB estimate of  error rate: 20%
## Confusion matrix:
##        forest rock snow class.error
## forest     63   14   13  0.30000000
## rock       15   52   27  0.44680851
## snow       16   15  285  0.09810127

## Predict the land cover based on random forest
landCover <- predict(snowMeadow, rf1)

## Plot output with map
par(mfrow=c(1,2))
plotRGB(snowMeadow)
plot(landCover, col=c("green","grey50","white"))

## Determine percent cover of snow
length(landCover[landCover==3])/sum( ## Percent snow cover = snow pixels / (snow pixels + bare ground pixels)
length(landCover[landCover==2]),
length(landCover[landCover==3]))

## [1] 0.6917295

Model validation for snow only

## create snow vs no snow
trainData[,"snowPA"] <- ifelse(trainData$class=="snow",1,0)
colnames(trainData)[3:5] <- c("Red","Green","Blue") ## revise names to RGB

 ## withhold 20% for sample testing
pres <- trainData[trainData$snowPA==1,]
abs <- trainData[trainData$snowPA==0,]
fold.p <- kfold(pres, k=5)
occtest.p <- pres[fold.p == 4, ]
occtrain.p <- pres[fold.p != 4, ]
fold.a <- kfold(abs, k=5)
occtest.a <-abs[fold.a == 4, ]
occtrain.a <- abs[fold.a != 4, ]
    
## Combine training
trainAll <- rbind(occtrain.p,occtrain.a)
testAll <- rbind(occtest.p,occtest.a)

## Evaluate the random forest
rf2 <- randomForest(snowPA ~ Red + Green + Blue, data=trainAll, importance=T)
rf2

## 
## Call:
##  randomForest(formula = snowPA ~ Red + Green + Blue, data = trainAll,      importance = T) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 1
## 
##           Mean of squared residuals: 0.1069126
##                     % Var explained: 54

varImpPlot(rf2)

eval1 <- evaluate(occtest.p, occtest.a, rf2)
eval1

## class          : ModelEvaluation 
## n presences    : 63 
## n absences     : 37 
## AUC            : 0.9086229 
## cor            : 0.7028173 
## max TPR+TNR at : 0.8067667

eval1@pcor

## [1] 3.625944e-16

## Determine threshold value
specSens <- threshold(eval1, stat="spec_sens")
specSens

## [1] 0.8067667

## confusion matrix
predictTest <- ifelse(predict(rf2, newdata =testAll)>specSens, 1,0)

conTable <- table(predicted = predictTest, observed = testAll$snowPA)
conTable

##          observed
## predicted  0  1
##         0 34 14
##         1  3 49

accuracy <- sum(diag(conTable))/nrow(testAll)
accuracy

## [1] 0.83

Compare estimated percentages with visual and snow pillow

## Rename raster
names(snow) <-  c("Red","Green","Blue") ## revise names to RGB

## set up Random forest
rf1 <- randomForest(class ~ Red + Green + Blue, data=trainAll, importance=T)
rf1

## 
## Call:
##  randomForest(formula = class ~ Red + Green + Blue, data = trainAll,      importance = T) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 1
## 
##         OOB estimate of  error rate: 21.5%
## Confusion matrix:
##        forest rock snow class.error
## forest     53   12    8   0.2739726
## rock       15   35   24   0.5270270
## snow       11   16  226   0.1067194

## Extract meadow cover for each meadow
getSnow <- function(j){
  meadowTemp <- crop(snow, j)  ## mask meadow from larger image
  meadowTemp <- mask(meadowTemp, j)
  meadowSnow <- predict(meadowTemp, rf1)
  ## Determine percent cover of snow
  snowPer <- freq(meadowSnow, value=3)/(freq(meadowSnow, value=3) + freq(meadowSnow, value=2)) ## snow / (snow+rock) i.e. omit forest
  return(snowPer)
}

## Snow percent
library(foreach)
library(doParallel)

cl <- makeCluster(6, type="PSOCK")
registerDoParallel(cores = cl)
meadowCover <- foreach(j = 1:17, .combine=c, .packages=c("raster","rgdal","randomForest")) %dopar% {
  getSnow(buffMeadow[buffMeadow$meadow == meadowList[j],])
}

## Snow cover
estSnow <- data.frame(year = 2012,  meadow = meadowList, snowCover=meadowCover)

## Compare with visual data

## Load visual data
Tempdata <- read.csv("data//snowData.csv", stringsAsFactors = F)
visual2012 <- Tempdata  %>% filter(year==2012) %>% select(meadow, percent) 
visual2012[visual2012$meadow=="g1","meadow"] <- "littleg"

## Join with modelled cover
bothEstimates <- merge(estSnow, visual2012, by="meadow")

## determine correlation
corVal <- cor(bothEstimates$snowCover, bothEstimates$percent)
corVal

## [1] 0.588996

## Plot patterns
ggplot(bothEstimates, aes(x=snowCover*100, y=percent)) + geom_point(size=3) + theme_Publication() + xlab("Random forest predicted snow cover (%)") + ylab("Visual estimate snow cover (%)") + annotate("text", x = 55, y = 95, label = "Meadows in 2012", size=8, hjust = 0) + annotate("text", x = 55, y = 90, label = paste0("cor = ",round(corVal,2)), size=6, hjust = 0)