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.
Ruta de imágenes
Mult<-brick("G:/Master_Script_R/Clasification/20190529_MULT_DoselZI.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/20190529_DSM_DoselZI.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/ZonaI.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 : 2497, 2928, 7311216, 7 (nrow, ncol, ncell, nlayers)
## resolution : 0.042, 0.042 (x, y)
## extent : 1078423, 1078546, 912214.1, 912319 (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.002349889, 0.011383231, 0.004699779, 0.020935378, 0.031555656, 0.124849314, 0.201388568
## max values : 0.08548104, 0.20086977, 0.12289616, 0.23408865, 0.41123065, 0.59011213, 0.84785229
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_ZI.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_ZI.shp", layer: "Training_Dosel_ZI"
## with 2 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)
## X20190529_MULT_DoselZI.1 X20190529_MULT_DoselZI.2 X20190529_MULT_DoselZI.3
## [1,] 0.01458762 0.05145342 0.01959259
## [2,] 0.01263447 0.05075151 0.01974517
## [3,] 0.02221714 0.06845197 0.02938888
## [4,] 0.03103685 0.06701762 0.03741512
## [5,] 0.02273594 0.06900129 0.02618448
## [6,] 0.01458762 0.04992752 0.01907378
## X20190529_MULT_DoselZI.4 X20190529_MULT_DoselZI.5
## [1,] 0.1406271 0.5157549
## [2,] 0.1372396 0.4651560
## [3,] 0.1621424 0.4657053
## [4,] 0.1493553 0.4564584
## [5,] 0.1555199 0.4799573
## [6,] 0.1289998 0.4399786
Tra_Sam <- rasterize(Tra_Pol, E_MULT); Tra_Sam
## class : RasterLayer
## dimensions : 2497, 2928, 7311216 (nrow, ncol, ncell)
## resolution : 0.042, 0.042 (x, y)
## extent : 1078423, 1078546, 912214.1, 912319 (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, 2 (min, max)
## attributes :
## ID LABEL Area
## 1 PE 170.327
## 2 PS 297.380
plot(Tra_Sam, axes=T, main="Distribucion de areas de entrenamiento")
names(Tra_Sam) <- c("LABEL")
Tra_Sam$LABEL
## class : RasterLayer
## dimensions : 2497, 2928, 7311216 (nrow, ncol, ncell)
## resolution : 0.042, 0.042 (x, y)
## extent : 1078423, 1078546, 912214.1, 912319 (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, 2 (min, max)
## attributes :
## ID LABEL Area
## 1 PE 170.327
## 2 PS 297.380
Class <- c("PE","PS")
ClassDF <- data.frame(ClassValue=c(1,2), ClassNames=Class)
# Hex codes of colors
ClassColor <- c("red","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
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 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 : 20000
## extent : 1078431, 1078546, 912215.9, 912319 (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 : 2795, 1
## max values : 7183234, 2
table(Samp$LABEL)
##
## 1 2
## 10000 10000
# Plot SpatialPointsDataFrame
#install.packages("rasterVis")
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
plot_SPDF <- levelplot(Tra_Sam$LABEL, col.regions=ClassColor,
main='Distribucion de areas de entrenamiento')
print(plot_SPDF + layer(sp.points(Samp, pch=3, cex=0.5, col=1)))
# Extract the layer values for the locations
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.02142367 0.06906233 0.03601129 0.17489891 0.4650645 263.9732
## 2 1 0.03140307 0.06543069 0.07608148 0.12463569 0.3122301 261.7375
## 3 1 0.01400778 0.03466850 0.01846342 0.09475853 0.4227970 265.1324
## 4 1 0.02124056 0.06594949 0.02938888 0.17660792 0.5130388 263.8622
## 5 1 0.01907378 0.06717022 0.02944991 0.17837797 0.4844739 263.4058
## 6 1 0.03418021 0.09228656 0.04550240 0.20263982 0.4959792 263.4763
## Chl_G Chl_RE Datt4 Maccioni MCARI OSAVI mND705
## 1 3.633230 2.034293 0.05498764 0.6762928 -0.0017073192 0.2857060 0.4859450
## 2 2.338386 1.176524 0.11669494 0.7943913 -0.0300324017 0.2746564 0.5015094
## 3 6.464348 3.571159 0.03014933 0.8113065 -0.0039158815 0.3371451 0.6700954
## 4 4.228598 2.347852 0.04589887 0.6956083 -0.0002603969 0.3031616 0.5198529
## 5 3.934121 2.189427 0.04555972 0.6727029 0.0004949529 0.2886646 0.4899853
## 6 2.785053 1.815398 0.06117395 0.6511754 0.0072414019 0.2786950 0.4654271
## NDRE REP_Li Vogelmann
## 1 0.4534096 712.9483 1.829524
## 2 0.4294097 708.2243 1.752571
## 3 0.6338227 707.5477 2.730918
## 4 0.4878308 712.1757 1.952480
## 5 0.4617864 713.0919 1.857998
## 6 0.4198847 713.9530 1.723795
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
#install.packages("rpart")
library(rpart)
## Warning: package 'rpart' was built under R version 4.0.5
# 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)
## Warning: package 'rpart.plot' was built under R version 4.0.5
rpart.plot(cart)
(PRED <- predict(E_All, cart, type='class'))
## class : RasterLayer
## dimensions : 2497, 2928, 7311216 (nrow, ncol, ncell)
## resolution : 0.042, 0.042 (x, y)
## extent : 1078423, 1078546, 912214.1, 912319 (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, 2 (min, max)
## attributes :
## ID value
## 1 1
## 2 2
(rat <- levels(PRED)[[1]])
## ID value
## 1 1 1
## 2 2 2
rat$Legn <- ClassDF$ClassNames; rat
## ID value Legn
## 1 1 1 PE
## 2 2 2 PS
levels(PRED) <- rat
levelplot(PRED, maxpixels = 1e6,
col.regions = ClassColor,
axes=FALSE)
EVALUACIÓN DE EXACTITUD DT
#install.packages("dismo")
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
set.seed(2021)
j <- kfold(SampData, k=2, by=SampData$ID_Class)
table(j)
## j
## 1 2
## 10000 10000
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
## 1 7400 2600
## 2 1294 8706
(cart_DT <- ovAcc(conmat))
## precision recall f1 OA kappa
## 1 0.8511617 0.7400 0.7916979 0.8053 0.6106
## 2 0.7700336 0.8706 0.8172346 0.8053 0.6106
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: 94%. Estimated remaining time: 1 seconds.
fit$variable.importance
## Blue Green Red RedEdge NIR DSM Chl_G Chl_RE
## 524.0822 545.2400 344.3681 296.7329 313.1013 419.0839 568.5443 409.8661
## Datt4 Maccioni MCARI OSAVI mND705 NDRE REP_Li Vogelmann
## 290.6318 871.7191 633.0671 1402.3349 1375.1807 582.3486 803.3722 619.8096
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)
## Predicting.. Progress: 92%. Estimated remaining time: 2 seconds.
## Predicting.. Progress: 80%. Estimated remaining time: 7 seconds.
## Predicting.. Progress: 84%. Estimated remaining time: 5 seconds.
## Predicting.. Progress: 67%. Estimated remaining time: 14 seconds.
classified
## class : RasterLayer
## dimensions : 2497, 2928, 7311216 (nrow, ncol, ncell)
## resolution : 0.042, 0.042 (x, y)
## extent : 1078423, 1078546, 912214.1, 912319 (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, 2 (min, max)
## attributes :
## ID value
## 1 1
## 2 2
plot(classified)
unique(classified$layer)
## [1] 1 2
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)
lc_colors <- c('red','green') %>%
colorRampPalette()
plot(classified, legend=TRUE, axes=TRUE, col=lc_colors(2))
region.brick <- crop(E_All, E)
region.lc <- crop(classified, E)
region.lc
## class : RasterLayer
## dimensions : 2497, 2928, 7311216 (nrow, ncol, ncell)
## resolution : 0.042, 0.042 (x, y)
## extent : 1078423, 1078546, 912214.1, 912319 (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, 2 (min, max)
## attributes :
## ID value
## 1 1
## 2 2
new_colors <- c('red','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(2))
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 : 20000
## extent : 1078431, 1078546, 912215.9, 912319 (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 : 2794, 1
## max values : 7186160, 2
table(TSamp$LABEL)
##
## 1 2
## 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 1 2 1 1 1 1 2 1 2 1 1 1 2 2 2 1 1 1 1 1 1 2 1 1 1 1 1 1 1
## Levels: 1 2
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
(conmat2 <- table(TSamp$prediction, TSamp$reference))
##
## 1 2
## 1 8588 1160
## 2 1412 8840
(cart_RF <- ovAcc(conmat2))
## precision recall f1 OA kappa
## 1 0.8588 0.8810012 0.8697590 0.8714 0.7428
## 2 0.8840 0.8622708 0.8730002 0.8714 0.7428