Clasificacion no Supervisada
Babotellos
31 de marzo de 2021
librerias
library(terra)
library(raster)entorno y procesamiento
input <- function(name,num='',type,foldata=F,folder){
ruta <- getwd()
if (foldata == T){
carpeta <- folder
file <- paste0(ruta,'/',carpeta,'/',name, num, type)
}else{
file <- paste0(ruta,'/',name, num, type)
}
}
bruto_landsat <- rast(input('LC08_L1TP_009057_20140820_20200911_02_T1_B',1:7,'.tif',T,'Datos'))
names(bruto_landsat) <- c('ultra-blue','blue','green','red','NIR','SWIR1','SWIR2')Recorte
xmx <- 374838
xmn <- 358682
ymx <- 467578
ymn <- 454847
e <- ext(xmn, xmx, ymn, ymx)
# crop landsat by the extent
landsat <- crop(bruto_landsat, e)
LSbrick <- brick(input('math',type = '.tif'))
landsat## class : SpatRaster
## dimensions : 424, 538, 7 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 358695, 374835, 454845, 467565 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : memory
## names : ultra-blue, blue, green, red, NIR, SWIR1, ...
## min values : 9496, 8530, 7553, 6579, 5577, 4680, ...
## max values : 28450, 31273, 33997, 36600, 48916, 65535, ...LSbrick## class : RasterBrick
## dimensions : 424, 538, 228112, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 358695, 374835, 454845, 467565 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : C:/Users/LENOVO/Desktop/11/RS/MyData/math.tif
## names : ultra.blue, blue, green, red, NIR, SWIR1, SWIR2
## min values : 9496, 8530, 7553, 6579, 5577, 4680, 4845
## max values : 28450, 31273, 33997, 36600, 48916, 65535, 65535no supervisadas
- Cuando no se tiene mucho conocimiento previo
Realizaremos una clasificacion sin supervision en un subconjunto espacial de la capa NDVI
ndvi <- ((LSbrick$NIR-LSbrick$red)/(LSbrick$NIR+LSbrick$red))
ndvi## class : RasterLayer
## dimensions : 424, 538, 228112 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 358695, 374835, 454845, 467565 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : memory
## names : layer
## values : -0.2194487, 0.5966747 (min, max)plot(ndvi,col=rev(terrain.colors(100)),main='landsat NDVI')
Clasificacion Kmeans
as matrix
nr <- values(ndvi)
str(nr) # concertir ndvi de spatraster a matrix## num [1:228112] 0.447 0.425 0.404 0.408 0.37 ...random locations and clusters
set.seed(99) #Random Locations
kmncluster <- kmeans(na.omit(nr),
centers = 10, #number clusters
iter.max = 500, #number of iterations
nstart = 5,
algorithm = 'Lloyd' #algorithm method
)
str(kmncluster)## List of 9
## $ cluster : int [1:228112] 2 6 6 6 10 10 10 5 5 10 ...
## $ centers : num [1:10, 1] 0.282 0.449 0.237 -0.043 0.325 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:10] "1" "2" "3" "4" ...
## .. ..$ : NULL
## $ totss : num 2523
## $ withinss : num [1:10] 4.71 4.29 4.62 3.12 4.69 ...
## $ tot.withinss: num 47.4
## $ betweenss : num 2475
## $ size : int [1:10] 29578 32524 25597 1659 30976 32355 18792 19155 8027 29449
## $ iter : int 212
## $ ifault : NULL
## - attr(*, "class")= chr "kmeans"matrix to SpatRaster
knr <- ndvi
values(knr) <- kmncluster$cluster
knr## class : RasterLayer
## dimensions : 424, 538, 228112 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 358695, 374835, 454845, 467565 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : memory
## names : layer
## values : 1, 10 (min, max)identify classes
mycolor <- c("greenyellow","green","forestgreen","darkblue",
"dodgerblue4","deepskyblue3","deepskyblue","purple4", "darkorchid1",
"magenta", "brown1")
plot(knr, main = 'Unsupervised classification', col = mycolor)
contrast data
#import vector data
valle <- shapefile(input('valle1','','.shp',T,'shapes'))#plot
vallecsv <- read.csv(input('TablaValle','','.csv',T,'shapes'))#plot
plot(knr, main = 'shapes', col = mycolor)
plot(valle, add=T)
muestreo uniforme parecido a rasterisar
library(sp)
samples <- spsample(valle, 10000, type='regular')
numero <- length(samples)
df_samples <- as(samples,"SpatialPointsDataFrame")
df_samples@data = data.frame(ID=1:numero, size=1)plot
plot(knr, main = 'sampling', col = mycolor)
plot(valle, add=T)
plot(df_samples, pch = 1, cex = (df_samples$size)/4, add=T)
extraer valores del raster en los puntos de muestreo
vallen <- valle[10:102,]
df_samples$estado <- over(df_samples, vallen)$estado
rasdf <- raster::extract(knr, df_samples)
head(rasdf)## [1] 1 3 3 1 3 8Valores por estado
library(modeest)
modams <- aggregate(rasdf, list(df_samples$estado),mfv)
listms <- aggregate(rasdf, list(df_samples$estado),list)
colnames(modams) <- c("estado","categoria")
colnames(listms) <- c("estado","categoria")
modamsclist <- listms$categoria[[1]]
montanaslist <- listms$categoria[[2]]
nlist <- listms$categoria[[3]]
plist <- listms$categoria[[4]]
riolist <- listms$categoria[[5]]
urbanolist <- listms$categoria[[6]]library(dplyr)
riocsv <- filter(vallecsv,estado == "rio")
urbanocsv <- filter(vallecsv,estado == "urbano")
desnudocsv <- filter(vallecsv,estado == "n")
pastocsv <- filter(vallecsv,estado == "p")
cultivocsv <- filter(vallecsv,estado == "c")
rioindex<- (riocsv[,1])
urbanoindex<- (urbanocsv[,1])
desnudoindex<- (desnudocsv[,1])
pastoindex<- (pastocsv[,1])
cultivoindex<- (cultivocsv[,1])Identificacion de cuerpos de agua
riotable<-aggregate(data.frame(Frecuencia=riolist),list(Categoria = riolist),length)
plot(riotable,type = "h", col="blue",lwd=4)
axis(side=1, at = seq(1,10,1), labels = seq(1,10,1))
colorrio <- c(rep('white',10))
colorrio[mfv(riolist)] <- "blue"
plot(knr, main = 'Water Bodies Using Unsupervised classification' , col = colorrio)
plot(valle[rioindex,],axes=T,xlim=c(358682,374838),ylim=c(454847,467578),add=T)
Identificacion de Cascos urbanos
urbanotable<-aggregate(data.frame(Frecuencia=urbanolist),
list(Categoria = urbanolist),length)
plot(urbanotable,type = "h", col="gray",lwd=4)
axis(side=1, at = seq(1,10,1), labels = seq(1,10,1))
colorurbano <- c(rep('white',10))
colorurbano[mfv(urbanolist)] <- "gray"
plot(knr, main = 'Urban Sites Using Unsupervised classification', col = colorurbano)
plot(valle[urbanoindex,],axes=T,xlim=c(358682,374838),ylim=c(454847,467578),add=T)
Identificacion de suelo desnudo
ntable<-aggregate(data.frame(Frecuencia=nlist),list(Categoria = nlist),length)
ntableo <- ntable[with(ntable, order(-ntable$Frecuencia)),]
plot(ntable,type = "h", col="brown",lwd=4)
axis(side=1, at = seq(1,10,1), labels = seq(1,10,1))
colordesnudo <- c(rep('white',10))
colordesnudo[ntableo$Categoria[1:4]] <- "brown"
plot(knr, main = 'Bare Ground Using Unsupervised classification', col = colordesnudo)
plot(valle[desnudoindex,],axes=T,xlim=c(358682,374838),ylim=c(454847,467578),add=T)
Identificacion de pasturas
ptable<-aggregate(data.frame(Frecuencia=plist),list(Categoria = plist),length)
ptableo <- ptable[with(ptable, order(-ptable$Frecuencia)),]
plot(ptable,type = "h", col="chartreuse",lwd=4)
axis(side=1, at = seq(1,10,1), labels = seq(1,10,1))
colorpasto <- c(rep('white',10))
colorpasto[ptableo$Categoria[1]] <- "chartreuse"
plot(knr, main = 'Grass Using Unsupervised classification', col = colorpasto)
plot(valle[pastoindex,],axes=T,xlim=c(358682,374838),ylim=c(454847,467578),add=T)
Identificacion de cultivo
ctable<-aggregate(data.frame(Frecuencia=clist),list(Categoria = clist),length)
plot(ctable,type = "h", col="darkgreen",lwd=4)
ctableo <- ctable[with(ctable, order(-ctable$Frecuencia)),]
axis(side=1, at = seq(1,10,1), labels = seq(1,10,1))
colorcultivo <- c(rep('white',10))
colorcultivo[ctableo$Categoria[1:3]] <- "darkgreen"
plot(knr, main = 'Crop Using Unsupervised classification', col = colorcultivo)
plot(valle[cultivoindex,],axes=T,xlim=c(358682,374838),ylim=c(454847,467578),add=T)
Clasificacion final
newcolor <- c(rep('white',10))
newcolor[mfv(riolist)] <- "blue"
newcolor[ntableo$Categoria[1:4]] <- "brown"
newcolor[ptableo$Categoria[1]] <- "chartreuse"
newcolor[ptableo$Categoria[2]] <- "greenyellow"
newcolor[ctableo$Categoria[1:2]] <- "darkgreen"
newcolor[ctableo$Categoria[3]] <- "chartreuse4"
plot(knr, main = 'Unsupervised classification cover', col = newcolor)
plot(valle[urbanoindex,],col='white', border='white',add=T)
plot(valle[98:99,],axes=T,xlim=c(358682,374838),ylim=c(454847,467578),col='white', border='white',add=T)
legend(x ='bottomright',legend=c('Water','Bare ground','low vegetation','grass','low crop','crop','undefined'), fill=c('blue','brown','greenyellow','chartreuse','chartreuse4','darkgreen','white'))
independiente de la clasificacion basada en los cluster del Kmeans la covertura va a ser la misma con este algoritmo
Kmean1
reproducibilidad
sessionInfo()## R version 4.0.4 (2021-02-15)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 17763)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=Spanish_Colombia.1252 LC_CTYPE=Spanish_Colombia.1252
## [3] LC_MONETARY=Spanish_Colombia.1252 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Colombia.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] dplyr_1.0.4 modeest_2.4.0 raster_3.3-7 sp_1.4-2 terra_1.0-10
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.5 pillar_1.4.6 compiler_4.0.4
## [4] highr_0.8 tools_4.0.4 timeSeries_3062.100
## [7] rpart_4.1-15 digest_0.6.25 tibble_3.0.2
## [10] jsonlite_1.7.2 evaluate_0.14 lifecycle_0.2.0
## [13] lattice_0.20-41 stable_1.1.4 clue_0.3-58
## [16] pkgconfig_2.0.3 rlang_0.4.10 DBI_1.1.0
## [19] yaml_2.2.1 rgdal_1.5-23 rmutil_1.1.5
## [22] xfun_0.21 statip_0.2.3 stringr_1.4.0
## [25] knitr_1.31 cluster_2.1.0 generics_0.1.0
## [28] vctrs_0.3.6 tidyselect_1.1.0 grid_4.0.4
## [31] glue_1.4.1 R6_2.4.1 spatial_7.3-13
## [34] rmarkdown_2.7 purrr_0.3.4 magrittr_1.5
## [37] ellipsis_0.3.1 fBasics_3042.89.1 codetools_0.2-18
## [40] htmltools_0.5.1.1 assertthat_0.2.1 stabledist_0.7-1
## [43] timeDate_3043.102 stringi_1.4.6 crayon_1.3.4