rs <- split(r)
trn <- st_extract(rs, pts)
trn$cls <- rep("land", 20)
trn$cls[c(8, 14, 15, 18, 19)] <- "water"
model <- MASS::lda(cls ~ ., st_drop_geometry(trn)) #linear discriminant analysis
pr <- predict(rs, model)
plot(pr[1], key.pos = 4, key.width = lcm(3.5), key.length = lcm(2))