library(raster)

## Loading required package: sp

library(randomForest)

## randomForest 4.6-14

## Type rfNews() to see new features/changes/bug fixes.

library(psych)

## 
## Attaching package: 'psych'

## The following object is masked from 'package:randomForest':
## 
##     outlier

setwd("~/Dropbox/BGU/1_Other/Amos_Bouskila/p_01_rgb_image_supervised_classification_example/")

###############################################

# Raster
r = brick("DSC_0588.JPG")

# Classified points
dat = read.csv("DSC_0588.csv", stringsAsFactors = FALSE)
dat$type = as.factor(dat$type)
xy = xyFromCell(r, dat$cell)

# Plot raster & classified points
plotRGB(r)
points(xy)
text(xy, as.character(dat$type), col = c("brown", "green", "yellow")[dat$type], pos = 3)

# Add 'greenness' index
dat$greenness = (dat$DSC_0588.2 - dat$DSC_0588.1) / (dat$DSC_0588.2 + dat$DSC_0588.1)

# Split to 'training' and 'test' data
test = sample(1:nrow(dat), 25)
test = (1:nrow(dat)) %in% test

# Fit model
fit = randomForest(type ~ DSC_0588.1 + DSC_0588.2 + DSC_0588.3 + greenness, data = dat[!test, ])
fit

## 
## Call:
##  randomForest(formula = type ~ DSC_0588.1 + DSC_0588.2 + DSC_0588.3 +      greenness, data = dat[!test, ]) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 2
## 
##         OOB estimate of  error rate: 18%
## Confusion matrix:
##       dry green soil class.error
## dry    14     5    1   0.3000000
## green   2    12    0   0.1428571
## soil    1     0   15   0.0625000

# Predict on 'test' data
dat$pred = predict(fit, dat)

# Evaluate model accuracy
tab = table(dat$type[test], dat$pred[test])
cohen.kappa(tab)

## Warning in cohen.kappa1(x, w = w, n.obs = n.obs, alpha = alpha, levels =
## levels): upper or lower confidence interval exceed abs(1) and set to +/- 1.

## Call: cohen.kappa1(x = x, w = w, n.obs = n.obs, alpha = alpha, levels = levels)
## 
## Cohen Kappa and Weighted Kappa correlation coefficients and confidence boundaries 
##                  lower estimate upper
## unweighted kappa  0.71     0.88     1
## weighted kappa    0.57     0.82     1
## 
##  Number of subjects = 25

# Predict on entire image
img = as.data.frame(r)
img$greenness = (img$DSC_0588.2 - img$DSC_0588.1) / (img$DSC_0588.2 + img$DSC_0588.1)
img$pred = predict(fit, img)
pred = r[[1]]
pred[] = as.numeric(factor(img$pred, levels = c("green", "dry", "soil")))

# Plot
plot(pred, col = c("brown", "green", "yellow"))