## ----images-------------------------------------------------------------------
workdir <- file.path(dirname(tempdir()))
datadir <- file.path(workdir, "reagro/sentinel")
dir.create(datadir, showWarnings=FALSE, recursive=TRUE)


## ---- message=FALSE-----------------------------------------------------------
ff <- agrodata::data_crop("sentinel", 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]
head(basename(fvv))


## ----spatrast, message=FALSE--------------------------------------------------
library(terra)
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 <- agrodata::data_crop("crop_ref")
head(cropref)
table(cropref$class)


## ----ref2---------------------------------------------------------------------
vecref <- vect(cropref, geom=c("x", "y"), crs=crs(vv))
vecref 
plot(vecref, "class", plg=list(x="bottomleft"))


## ----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[[1:4]])


## ----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[,-1])
ref$class <- as.factor(ref$class)
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-----------------------------------------
levels(rf.pred) <- data.frame(ID=1:4, class=c("Grassland", "Maize", "Pasture", "Wheat"))
plot(rf.pred, col=c("green", "yellow", "darkgreen", "magenta"))


## ----sample-------------------------------------------------------------------
library(predicts)
set.seed(530)
k <- sample(5, nrow(ref), replace=TRUE)
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------------------------------------------------------------------
library(luna)
# overall accuracy
evaluate(conmat, "overall")
# kappa
evaluate(conmat, "kappa")
# user and producer accuracy
evaluate(conmat, "class")