Librerías
library(raster)
## Loading required package: sp
library(sp)
library(rgdal)
## rgdal: version: 1.5-19, (SVN revision 1092)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: C:/Users/user/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: C:/Users/user/Documents/R/win-library/4.0/rgdal/proj
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
library(rasterVis)
## Warning: package 'rasterVis' was built under R version 4.0.5
## Loading required package: terra
## Warning: package 'terra' was built under R version 4.0.4
## terra version 1.0.10
##
## Attaching package: 'terra'
## The following object is masked from 'package:rgdal':
##
## project
## Loading required package: lattice
## Loading required package: latticeExtra
## Warning: package 'latticeExtra' was built under R version 4.0.4
library(rpart)
## Warning: package 'rpart' was built under R version 4.0.5
library(rpart.plot)
## Warning: package 'rpart.plot' was built under R version 4.0.5
library(dismo)
## Warning: package 'dismo' was built under R version 4.0.5
##
## Attaching package: 'dismo'
## The following object is masked from 'package:terra':
##
## voronoi
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.0.4
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:terra':
##
## collapse, desc, intersect, near, select, union
## The following objects are masked from 'package:raster':
##
## intersect, select, union
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(RColorBrewer)
Ruta de imágenes
Mult<-brick("G:/Master_Script_R/Clasification/20190619_MULT_DoselZII.tif")
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition
DSM<-raster("G:/Master_Script_R/Clasification/20190619_DSM_DoselZII.tif")
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition
E <- shapefile("G:/Master_Script_R/Clasification/ZonaII.shp")
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS
## = dumpSRS, : Discarded datum Marco_Geocentrico_Nacional_de_Referencia in CRS
## definition: +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1
## +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
par(mfrow = c(1,1))
plotRGB(Mult, r=3,g=2,b=1, axes=T, stretch="lin")
plot(E, lwd=3, add=TRUE)
E_MULT <- crop(Mult, E);
E_DSM <- crop(DSM, E);
E_PLC<-stack(E_MULT,E_DSM)
names(E_PLC) <- c('Blue','Green','Red','RedEdge','NIR','DSM')
Índices multiespectrales
rRE=(E_PLC[[3]]+E_PLC[[4]])/2
rNIR=(E_PLC[[4]]+E_PLC[[5]])/2
ORTHOM7<-stack(E_PLC[[1]],E_PLC[[2]],E_PLC[[3]],rRE,
E_PLC[[4]],rNIR,E_PLC[[5]])
# Change band names
names(ORTHOM7) <- c('Blue','Green','Red','rRE','RedEdge','rNIR','NIR')
ORTHOM7
## class : RasterStack
## dimensions : 2920, 2917, 8517640, 7 (nrow, ncol, ncell, nlayers)
## resolution : 0.042, 0.042 (x, y)
## extent : 1078423, 1078546, 912196.1, 912318.7 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs
## names : Blue, Green, Red, rRE, RedEdge, rNIR, NIR
## min values : 0.003906310, 0.012298772, 0.004943923, 0.025589379, 0.043244068, 0.116166936, 0.186861977
## max values : 0.06863508, 0.15603876, 0.12857252, 0.18913558, 0.32947281, 0.49941254, 0.69748992
CARTER4 <- function(IMG, x,y) {
RE<-IMG[[x]]
NIR<-IMG[[y]]
CARTER4<- (RE/NIR)
return(CARTER4)
}
Chl_G <- function(IMG, x,y) {
G<-IMG[[x]]
rNIR<-IMG[[y]]
Chl_G<- (rNIR/G)-1
return(Chl_G)
}
Chl_RE <- function(IMG, x,y) {
rRE<-IMG[[x]]
rNIR<-IMG[[y]]
Chl_RE<- (rNIR/rRE)-1
return(Chl_RE)
}
Datt <- function(IMG, x,y,z) {
R<-IMG[[x]]
RE<-IMG[[y]]
NIR<-IMG[[z]]
Datt<- (NIR-RE)/(NIR-R)
return(Datt)
}
Datt2 <- function(IMG, x,y) {
RE<-IMG[[x]]
NIR<-IMG[[y]]
Datt2<- (NIR/RE)
return(Datt2)
}
Datt4 <- function(IMG, x,y,z) {
G<-IMG[[x]]
R<-IMG[[y]]
rRE<-IMG[[z]]
Datt4<- R/G*rRE
return(Datt4)
}
Maccioni <- function(IMG, x,y,z) {
rRE<-IMG[[x]]
RE<-IMG[[y]]
rNIR<-IMG[[z]]
Maccioni<-(rNIR-RE)/(rNIR-rRE)
return(Maccioni)
}
MCARI <- function(IMG, x,y,z,v,w) {
G<-IMG[[x]]
R<-IMG[[y]]
rRE<-IMG[[z]]
RE<-IMG[[v]]
rNIR<-IMG[[w]]
MCARI<- ((rRE-R)-2*(rRE-G))*(rNIR-RE)/(rNIR-R)
return(MCARI)
}
OSAVI <- function(IMG, x,y,z,v) {
R<-IMG[[x]]
RE<-IMG[[y]]
rNIR<-IMG[[z]]
NIR<-IMG[[v]]
OSAVI<-(((1+0.16)*(rNIR-RE))/(NIR-R+0.16))
return(OSAVI)
}
mND705 <- function(IMG, x,y,z) {
B<-IMG[[x]]
RE<-IMG[[y]]
NIR<-IMG[[z]]
mND705<-(NIR-RE)/(NIR+RE-2*B)
return(mND705)
}
NDRE <- function(IMG, x,y) {
RE<-IMG[[x]]
NIR<-IMG[[y]]
NDRE<- (NIR-RE)/(NIR+RE)
return(NDRE)
}
REP_Li <- function(IMG, x,y,z) {
rRE<-IMG[[x]]
RE<-IMG[[y]]
rNIR<-IMG[[z]]
REP_Li<-700+40*(RE-rRE)/(rNIR-rRE)
return(REP_Li)
}
Vogelmann <- function(IMG, x,y) {
RE<-IMG[[x]]
rNIR<-IMG[[y]]
Vogelmann<- (rNIR/RE)
return(Vogelmann)
}
VogelmannX <- function(IMG, x,y,z) {
RE<-IMG[[x]]
rNIR<-IMG[[y]]
NIR<-IMG[[z]]
VogelmannX<-(rNIR-NIR)/(RE+rNIR)
return(VogelmannX)
}
MS_CARTER4=CARTER4(ORTHOM7,5,7)
MS_Chl_G=Chl_G(ORTHOM7,2,6)
MS_Chl_RE=Chl_RE(ORTHOM7,4,6)
MS_Datt=Datt(ORTHOM7,3,5,7)
MS_Datt2=Datt2(ORTHOM7,5,7)
MS_Datt4=Datt4(ORTHOM7,2,3,4)
MS_Maccioni=Maccioni(ORTHOM7,4,5,6)
MS_MCARI=MCARI(ORTHOM7,2,3,4,5,6)
MS_OSAVI=OSAVI(ORTHOM7,3,5,6,7)
MS_MCARI.OSAVI = MS_MCARI/MS_OSAVI
MS_mND705=mND705(ORTHOM7,1,5,7)
MS_NDRE=NDRE(ORTHOM7,5,7)
MS_REP_Li=REP_Li(ORTHOM7,4,5,6)
MS_Vogelmann=Vogelmann(ORTHOM7,5,6)
MS_VogelmannX=VogelmannX(ORTHOM7,5,6,7)
MS_Index<-stack(MS_Chl_G,MS_Chl_RE,MS_Datt4,MS_Maccioni,
MS_MCARI,MS_OSAVI,MS_mND705,MS_NDRE,
MS_REP_Li,MS_Vogelmann)
names(MS_Index) <- c('Chl_G','Chl_RE','Datt4','Maccioni',
'MCARI','OSAVI','mND705','NDRE',
'REP_Li','Vogelmann')
plot(E_PLC, col=gray(0:100/100))
plot(MS_Index, col=gray(0:100/100))
E_All <- stack(E_PLC,MS_Index)
Áreas de entrenamiento
par(mfrow=c(1,1))
Tra_Pol<-readOGR('G:/Master_Script_R/Clasification/Training_Dosel_ZII.shp')
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS
## = dumpSRS, : Discarded datum Marco_Geocentrico_Nacional_de_Referencia in CRS
## definition: +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1
## +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## OGR data source with driver: ESRI Shapefile
## Source: "G:\Master_Script_R\Clasification\Training_Dosel_ZII.shp", layer: "Training_Dosel_ZII"
## with 3 features
## It has 2 fields
plotRGB(E_MULT, r=3,g=2,b=1, axes=T, stretch="lin")
plot(Tra_Pol, axes=T, lwd=3, add=TRUE)
# Generate 10000 point samples from the polygons
Pt_Samp <- spsample(Tra_Pol, 10000, type='regular')
## Warning in proj4string(obj): CRS object has comment, which is lost in output
# Add the land cover class to the points
Pt_Samp$LABEL <- over(Pt_Samp, Tra_Pol)$LABEL
# extract values with points
DF <- extract(E_MULT, Pt_Samp)
# To see some of the reflectance values
head(DF)
## X20190619_MULT_DoselZII.1 X20190619_MULT_DoselZII.2
## [1,] 0.01629664 0.04303044
## [2,] 0.02108797 0.04727245
## [3,] 0.01882963 0.04501411
## [4,] NA NA
## [5,] NA NA
## [6,] 0.01886015 0.04724193
## X20190619_MULT_DoselZII.3 X20190619_MULT_DoselZII.4
## [1,] 0.01995880 0.1183185
## [2,] 0.02212558 0.1205768
## [3,] 0.01986725 0.1221942
## [4,] NA NA
## [5,] NA NA
## [6,] 0.02124056 0.1368429
## X20190619_MULT_DoselZII.5
## [1,] 0.3783627
## [2,] 0.3834897
## [3,] 0.3801633
## [4,] NA
## [5,] NA
## [6,] 0.4201724
Tra_Sam <- rasterize(Tra_Pol, E_MULT); Tra_Sam
## class : RasterLayer
## dimensions : 2920, 2917, 8517640 (nrow, ncol, ncell)
## resolution : 0.042, 0.042 (x, y)
## extent : 1078423, 1078546, 912196.1, 912318.7 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## source : memory
## names : layer
## values : 1, 3 (min, max)
## attributes :
## ID LABEL Area
## 1 PE 107.0300
## 2 PI 44.7418
## 3 PS 216.8120
plot(Tra_Sam, axes=T, main="Distribucion de areas de entrenamiento")
names(Tra_Sam) <- c("LABEL")
Tra_Sam$LABEL
## class : RasterLayer
## dimensions : 2920, 2917, 8517640 (nrow, ncol, ncell)
## resolution : 0.042, 0.042 (x, y)
## extent : 1078423, 1078546, 912196.1, 912318.7 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## source : memory
## names : LABEL
## values : 1, 3 (min, max)
## attributes :
## ID LABEL Area
## 1 PE 107.0300
## 2 PI 44.7418
## 3 PS 216.8120
Class <- c("PE","PI","PS")
ClassDF <- data.frame(ClassValue=c(1,2,3), ClassNames=Class)
# Hex codes of colors
ClassColor <- c("red","gray","green")
# Now we ratify (RAT="Raster Attribute Table") the Class (define RasterLayer as a categorical variable)
LABEL <- Tra_Sam[[1]]
LABEL <- ratify(LABEL)
RAT <- levels(LABEL)[[1]]; RAT
## ID
## 1 1
## 2 2
## 3 3
RAT$ID <- Class
levels(LABEL) <- RAT; RAT
## Warning in .checkLevels(levels(x)[[1]], value): the values in the "ID" column in
## the raster attributes (factors) data.frame have changed
## Warning in .checkLevels(levels(x)[[1]], value): NAs introducidos por coerción
## ID
## 1 PE
## 2 PI
## 3 PS
plot(Tra_Sam, axes=T, col=ClassColor,
main='Distribucion de areas de entrenamiento')
# Set the random number generator to reproduce the results
set.seed(100)
# Sampling
Samp <- sampleStratified(LABEL, size=10000, na.rm=TRUE, sp=TRUE)
Samp
## class : SpatialPointsDataFrame
## features : 30000
## extent : 1078425, 1078545, 912215.3, 912318.6 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## variables : 2
## names : cell, LABEL
## min values : 5699, 1
## max values : 7181768, 3
table(Samp$LABEL)
##
## 1 2 3
## 10000 10000 10000
# Plot SpatialPointsDataFrame
library(rasterVis)
plot_SPDF <- levelplot(Tra_Sam$LABEL, col.regions=ClassColor,
main='Distribucion of Training Sites')
print(plot_SPDF + layer(sp.points(Samp, pch=3, cex=0.5, col=1)))
# Extract the layer values for the locations
# To keep the spatial information you use `sp=TRUE` argument in the `extract` function
SampVals <- extract(E_All, Samp, df=TRUE)
SampVals <- SampVals[,-1] # drop the ID column
# Combine the class information with extracted values
SampData <- data.frame(ID_Class = Samp$LABEL, SampVals)
SampData$ID_Class <- as.factor(SampData$ID_Class)
head(SampData)
## ID_Class Blue Green Red RedEdge NIR DSM
## 1 1 0.04037537 0.11511406 0.04452582 0.2700236 0.5219806 268.0040
## 2 1 0.04126040 0.07998779 0.05502403 0.1890898 0.4782178 268.7221
## 3 1 0.02124056 0.09094377 0.02688640 0.2256504 0.4944533 268.4017
## 4 1 0.02468910 0.09057756 0.02801557 0.2647440 0.5706264 267.4648
## 5 1 0.03384451 0.12832837 0.03427176 0.3159533 0.5990997 267.5320
## 6 1 0.03863584 0.09213398 0.03622492 0.2417029 0.5609522 267.5410
## Chl_G Chl_RE Datt4 Maccioni MCARI OSAVI mND705
## 1 2.440085 1.517901 0.06083346 0.5277086 0.010189199 0.2292477 0.3542436
## 2 3.171309 1.733592 0.08396361 0.6832047 -0.008874928 0.2875446 0.4944160
## 3 2.959060 1.851480 0.03732969 0.5748972 0.011591002 0.2484288 0.3966853
## 4 3.611354 1.853435 0.04527515 0.5637233 0.002653128 0.2525036 0.3891671
## 5 2.565280 1.612757 0.04676608 0.5012967 0.015812013 0.2265709 0.3341497
## 6 3.355913 1.887998 0.05463735 0.6084099 0.003969436 0.2704209 0.4401111
## NDRE REP_Li Vogelmann
## 1 0.3181258 718.8917 1.466546
## 2 0.4332754 712.6718 1.764525
## 3 0.3732836 717.0041 1.595618
## 4 0.3661637 717.4511 1.577694
## 5 0.3094317 719.9481 1.448083
## 6 0.3977415 715.6636 1.660417
SampData.df <- na.omit(SampData)
(predictors <- colnames(SampData.df)[2:17])
## [1] "Blue" "Green" "Red" "RedEdge" "NIR" "DSM"
## [7] "Chl_G" "Chl_RE" "Datt4" "Maccioni" "MCARI" "OSAVI"
## [13] "mND705" "NDRE" "REP_Li" "Vogelmann"
(fo <- as.formula(paste("ID_Class ~", paste(predictors, collapse = "+"))))
## ID_Class ~ Blue + Green + Red + RedEdge + NIR + DSM + Chl_G +
## Chl_RE + Datt4 + Maccioni + MCARI + OSAVI + mND705 + NDRE +
## REP_Li + Vogelmann
ovAcc <- function(conmat)
{
# number of cases
n = sum(conmat)
# number of correctly classified cases per class
diag = diag(conmat)
# Overall Accuracy
OA = sum(diag) / n
# observed (true) cases per class
rowsums = apply(conmat, 1, sum)
p = rowsums / n
# predicted cases per class
colsums = apply(conmat, 2, sum)
q <- colsums / n
expAccuracy <- sum(p*q)
kappa <- (OA - expAccuracy) / (1 - expAccuracy)
# Producer accuracy
PA <- diag / colsums
# User accuracy
UA <- diag / rowsums
outAcc <- data.frame(producerAccuracy = PA, userAccuracy = UA)
# Overall Accuracy & KAPPA
precision = diag / colsums
recall = diag / rowsums
f1 = 2 * precision * recall / (precision + recall)
(global_acc = data.frame(precision, recall, f1, OA, kappa))
}
Decision Tree
#library(rpart)
# Train the model
cart <- rpart(as.factor(ID_Class)~., data=SampData, method='class', minsplit=5)
# print(model.class)
# Plot the trained classification tree
#library(rpart.plot)
rpart.plot(cart)
(PRED <- predict(E_All, cart, type='class'))
## class : RasterLayer
## dimensions : 2920, 2917, 8517640 (nrow, ncol, ncell)
## resolution : 0.042, 0.042 (x, y)
## extent : 1078423, 1078546, 912196.1, 912318.7 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs
## source : memory
## names : layer
## values : 1, 3 (min, max)
## attributes :
## ID value
## 1 1
## 2 2
## 3 3
(rat <- levels(PRED)[[1]])
## ID value
## 1 1 1
## 2 2 2
## 3 3 3
rat$Legn <- ClassDF$ClassNames; rat
## ID value Legn
## 1 1 1 PE
## 2 2 2 PI
## 3 3 3 PS
levels(PRED) <- rat
levelplot(PRED, maxpixels = 1e6,
col.regions = ClassColor,
axes=FALSE)
EVALUACIÓN DE EXACTITUD DT
#library(dismo)
set.seed(2021)
j <- kfold(SampData, k=3, by=SampData$ID_Class)
table(j)
## j
## 1 2 3
## 9999 10002 9999
x <- list()
for (k in 1:2) {
train <- SampData[j!= k, ]
test <- SampData[j == k, ]
cart <- rpart(as.factor(ID_Class)~., data=train, method='class', minsplit=5)
pclass <- predict(cart, test, type='class')
# create a data.frame using the reference and prediction
x[[k]] <- cbind(test$ID_Class, as.integer(pclass))
}
y <- do.call(rbind, x)
y <- data.frame(y)
colnames(y) <- c('observed', 'predicted')
(conmat <- table(y))
## predicted
## observed 1 2 3
## 1 5927 51 689
## 2 712 5951 4
## 3 789 13 5865
(cart_DT <- ovAcc(conmat))
## precision recall f1 OA kappa
## 1 0.7979268 0.8890055 0.8410074 0.8871056 0.8306585
## 2 0.9893599 0.8926054 0.9384955 0.8871056 0.8306585
## 3 0.8943275 0.8797060 0.8869565 0.8871056 0.8306585
Random Forest
library(ranger)
## Warning: package 'ranger' was built under R version 4.0.5
fit <- ranger(fo, data = SampData.df, importance="impurity")
## Growing trees.. Progress: 77%. Estimated remaining time: 9 seconds.
fit$variable.importance
## Blue Green Red RedEdge NIR DSM Chl_G Chl_RE
## 432.9994 628.5973 1487.0434 416.2064 1592.1743 751.0039 537.1216 1263.2825
## Datt4 Maccioni MCARI OSAVI mND705 NDRE REP_Li Vogelmann
## 2763.2239 816.4319 3364.3300 2316.9535 1148.4591 733.7941 966.5803 781.1498
fo
## ID_Class ~ Blue + Green + Red + RedEdge + NIR + DSM + Chl_G +
## Chl_RE + Datt4 + Maccioni + MCARI + OSAVI + mND705 + NDRE +
## REP_Li + Vogelmann
classified <- predict(E_All, fit,
fun=function(model, ...) predict(model, ...)$predictions)
classified
## class : RasterLayer
## dimensions : 2920, 2917, 8517640 (nrow, ncol, ncell)
## resolution : 0.042, 0.042 (x, y)
## extent : 1078423, 1078546, 912196.1, 912318.7 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs
## source : memory
## names : layer
## values : 1, 3 (min, max)
## attributes :
## ID value
## 1 1
## 2 2
## 3 3
plot(classified)
unique(classified$layer)
## [1] 1 2 3
library(dplyr)
library(RColorBrewer)
lc_colors <- c('red','gray','green') %>%
colorRampPalette()
plot(classified, legend=TRUE, axes=TRUE, col=lc_colors(3))
region.brick <- crop(E_All, E)
region.lc <- crop(classified, E)
region.lc
## class : RasterLayer
## dimensions : 2920, 2917, 8517640 (nrow, ncol, ncell)
## resolution : 0.042, 0.042 (x, y)
## extent : 1078423, 1078546, 912196.1, 912318.7 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +units=m +no_defs
## source : memory
## names : layer
## values : 1, 3 (min, max)
## attributes :
## ID value
## 1 1
## 2 2
## 3 3
new_colors <- c('red','gray','green') %>%
colorRampPalette()
par(mfrow=c(1,2))
plotRGB(region.brick, r = 3, g = 2, b = 1, stretch = "lin", axes=T)
plot(region.lc, legend=F, axes=F, col=lc_colors(3))
EVALUACIÓN DE EXACTITUD RF
set.seed(2021)
# Sampling
TSamp <- sampleStratified(Tra_Sam$LABEL, size=10000, na.rm=TRUE, sp=TRUE)
TSamp
## class : SpatialPointsDataFrame
## features : 30000
## extent : 1078425, 1078545, 912215.3, 912318.6 (xmin, xmax, ymin, ymax)
## crs : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## variables : 2
## names : cell, LABEL
## min values : 5692, 1
## max values : 7181770, 3
table(TSamp$LABEL)
##
## 1 2 3
## 10000 10000 10000
TSamp$reference <- as.factor(TSamp$LABEL)
prediction <- raster::extract(classified, TSamp)
TSamp$prediction <- as.factor(prediction)
TSamp$prediction[700:730]
## [1] 1 1 3 1 1 1 1 1 3 1 1 1 1 1 1 3 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1
## Levels: 1 2 3
TSamp$reference[700:730]
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## Levels: 1 2 3
(conmat2 <- table(TSamp$prediction, TSamp$reference))
##
## 1 2 3
## 1 9222 126 564
## 2 103 9872 5
## 3 675 2 9431
(cart_RF <- ovAcc(conmat2))
## precision recall f1 OA kappa
## 1 0.9222 0.9303874 0.9262756 0.9508333 0.92625
## 2 0.9872 0.9891784 0.9881882 0.9508333 0.92625
## 3 0.9431 0.9330233 0.9380346 0.9508333 0.92625