## ---- global_options, echo=FALSE, message=FALSE--------------------------
library(knitr)
library(agrin)
library(randomForest)


## ----images--------------------------------------------------------------
workdir <- file.path(dirname(tempdir()))
datadir <- file.path(workdir, "agrin/sentinel")
dir.create(datadir, showWarnings = FALSE, recursive=TRUE)
library(agrin)
datadir
ff <- crop_data("sentinel1", datadir)
head(ff)


## ----haveimg-------------------------------------------------------------
imgfile <- file.path(workdir, "img.tif")
imgexists <- file.exists(imgfile)


## ----polvv---------------------------------------------------------------
fvv <- grep("VV", ff, value = TRUE)
i <- order(substr(basename(fvv), 15, 27))
fvv <- fvv[i]
fvv

vv <- rast(fvv)
# make prettier (shorter) layer names
names(vv) <- substr(names(vv), 15, 27)
names(vv)
vv


## ----polvh---------------------------------------------------------------
fvh <- grep("VH", ff, value = TRUE)
i <- order(substr(basename(fvh), 15, 27))
fvh <- fvh[i]
vh <- rast(fvh)
names(vh) <- substr(names(vh), 15, 27)
names(vh)


## ----ref1----------------------------------------------------------------
cropref <- crop_data("crop_ref")
head(cropref)
table(cropref$class)


## ----ref2----------------------------------------------------------------
vecref <- vect(cropref[,1:2], att=cropref, crs=crs(vv))
vecref 
spplot(vecref, "class")


## ----features------------------------------------------------------------
compFeatures <- function(x, name){
	# standard deviation
	stdev <- stdev(x)
	# quantiles
	quant <- app(x, function(i) quantile(i, c(.1, .25, .5, .75, .9)))
	# trend
	n <- nlyr(x)
	trend <- 100 * (mean(x[[1:5]] - mean(x[[(n-4):n]])))

	feat <- c(quant, stdev, trend)
	names(feat) <- paste0(name, c(".1", ".25", ".5", ".75", ".9", ".sd", "trend"))
	return(feat)
} 


## ----vvfeatures, fig.width=12, fig.height=12-----------------------------
fvv <- file.path(workdir, "vv_ft.tif")

if (!file.exists(fvv)) {
	vv_ft <- compFeatures(vv, "vv")
	writeRaster(vv_ft, fvv)
} else {
	vv_ft <- rast(fvv)
}

vv_ft
plot(vv_ft)


## ----vhfeatures, fig.width=12, fig.height=8------------------------------
fvh <- file.path(workdir, "vh_ft.tif")
if (!file.exists(fvh)) {
	sd <- stdev(vh)
	mn <- mean(vh)
	vh_ft <- c(mn, sd)
	names(vh_ft) <- c("mean", "stdev")
	writeRaster(vh_ft, fvh)
} else {
	vh_ft <-- rast(fvh)
}
plot(vh_ft)


## ----scaling-------------------------------------------------------------
img <- c(vv_ft, vh_ft)


## ----getref--------------------------------------------------------------
ref <- extract(img, vecref, drop=TRUE)
ref <- data.frame(class=cropref$class, ref)
head(ref)


## ----RFtrain-------------------------------------------------------------
library(randomForest)
set.seed(888)
rf.m <- randomForest(class~., data=ref, ntree=250, importance=TRUE)


## ----RFprint-------------------------------------------------------------
rf.m


## ----RFplot--------------------------------------------------------------
cols <- c("green", "yellow", "darkgreen", "magenta", "blue")
plot(rf.m, col=cols, lwd=3)
legend("topright", col=cols , c("OOB", "Grass", "Maize","Pasture","Wheat"), lty=1, lwd=3)


## ----RFPlots-------------------------------------------------------------
varImpPlot(rf.m, main="")


## ----RFplotsdata---------------------------------------------------------
importance(rf.m)


## ----partialplot---------------------------------------------------------
partialPlot(rf.m, ref, "vvtrend", "Maize")


## ----RFpredict-----------------------------------------------------------
rf.pred <- predict(img, rf.m)
rf.pred


## ----RFmap, fig.width=6, fig.height=6------------------------------------
rfp <- as.factor(rf.pred)
levels(rfp) <- c("Grassland", "Maize", "Pasture", "Wheat")
plot(rfp, col=c("green", "yellow", "darkgreen", "magenta"))


## ----sample--------------------------------------------------------------
library(luna)
set.seed(530)
k <- kfold(ref, 5)
nrow(ref) / 5
table(k)


## ----split---------------------------------------------------------------
train <- ref[k != 1, ]
test  <- ref[k == 1, ]


## ----predfit-------------------------------------------------------------
rf <- randomForest(class~., data=train, ntree=250, importance=TRUE)


## ----predtest------------------------------------------------------------
pred <- predict(rf, test)


## ----accuracy------------------------------------------------------------
conmat <- table(observed=test$class, predicted=pred)
conmat


## ----corclas-------------------------------------------------------------
# overall accuracy
evaluate(conmat, "overall")
# kappa
evaluate(conmat, "kappa")
# user and producer accuracy
evaluate(conmat, "class")