En este cuaderno se realiza una clasificación no supervisada (UC), el cual se caracteriza por la agrupación de pixeles que poseen regiones del espectro visible que pueden asociarse con diferentes clases de coberturas del suelo. La clasificación no supervisada K-means agrupa pixeles con misma respuesta espectral posteriormente idéntifica el centroide y su valor para cada grupo para luego dividir el espacio entre las celdas, es decir, cada observación como una única agrupación.

El ejercicio imagen satelital que pertenece al estado de Illionis en estados unidos: Landsat 8 "LC08_L1TP_022031_20200831_20200906_01_T1.TIF

Se carga la ruta de trabajo y se visualiza:

setwd('C:\\Users\\Antonio\\Desktop\\PERCEPCION_REMOTA\\TALLER\\LABORATORIO_5')
getwd()
## [1] "C:/Users/Antonio/Desktop/PERCEPCION_REMOTA/TALLER/LABORATORIO_5"

Se verifica la lista de archivos que están en la ruta:

list.files()
##  [1] "Lab-5-Supervised-classification.html"           
##  [2] "Lab-5-Supervised-classification.Rmd"            
##  [3] "Lab 5 Supervised classification.Rmd"            
##  [4] "LC08_L1TP_022031_20200831_20200906_01_T1_B1.TIF"
##  [5] "LC08_L1TP_022031_20200831_20200906_01_T1_B2.TIF"
##  [6] "LC08_L1TP_022031_20200831_20200906_01_T1_B3.TIF"
##  [7] "LC08_L1TP_022031_20200831_20200906_01_T1_B4.TIF"
##  [8] "LC08_L1TP_022031_20200831_20200906_01_T1_B5.TIF"
##  [9] "LC08_L1TP_022031_20200831_20200906_01_T1_B6.TIF"
## [10] "LC08_L1TP_022031_20200831_20200906_01_T1_B7.TIF"

Se usa la función paste0 para formar el nombre de la banda 2 a la banda 7.

rfiles <- paste0('LC08_L1TP_022031_20200831_20200906_01_T1_B',2:7,".TIF")
rfiles
## [1] "LC08_L1TP_022031_20200831_20200906_01_T1_B2.TIF"
## [2] "LC08_L1TP_022031_20200831_20200906_01_T1_B3.TIF"
## [3] "LC08_L1TP_022031_20200831_20200906_01_T1_B4.TIF"
## [4] "LC08_L1TP_022031_20200831_20200906_01_T1_B5.TIF"
## [5] "LC08_L1TP_022031_20200831_20200906_01_T1_B6.TIF"
## [6] "LC08_L1TP_022031_20200831_20200906_01_T1_B7.TIF"

Con los archivos “rfiles” se se crea un nuevo raster rast y se guardan en la variable landsat

landsat <- rast(rfiles)
names(landsat) <- c('blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')
landsat
## class       : SpatRaster 
## dimensions  : 7871, 7751, 6  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 403485, 636015, 4504485, 4740615  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs 
## sources     : LC08_L1TP_022031_20200831_20200906_01_T1_B2.TIF  
##               LC08_L1TP_022031_20200831_20200906_01_T1_B3.TIF  
##               LC08_L1TP_022031_20200831_20200906_01_T1_B4.TIF  
##               ... and 3 more source(s)
## names       : blue, green, red, NIR, SWIR1, SWIR2
plotRGB(landsat, r=6, g=5, b=3,  axes=TRUE, colNA="white", stretch='lin', main="Landsat 8 Falso color (para zonas urbanas)")

Calculo del NDVI (Índice de Vegetación de Diferencia Normalizada)

ndvi <- (landsat[['NIR']] - landsat[['red']]) / (landsat[['NIR']] + landsat[['red']])
ndvi
## class       : SpatRaster 
## dimensions  : 7871, 7751, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 403485, 636015, 4504485, 4740615  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs 
## source      : spat_XhQmHSDDSvNOrWD.tif 
## name        : NIR 
## min value   :  -1 
## max value   :   1

A partir del cálculo de NDVI se empleara el método de clasificación no supervisión “KMEANS” en el que se agruparan los datos.

Se carga un shpaefile para realizar un corte con el fin de tener una rea de estudio especifica:

knitr::opts_chunk$set(echo = TRUE)
setwd("C:\\Users\\Antonio\\Desktop\\PERCEPCION_REMOTA\\TALLER\\LABORATORIO 4\\shapefile_corte")
getwd()
## [1] "C:/Users/Antonio/Desktop/PERCEPCION_REMOTA/TALLER/LABORATORIO_5"
shape<-("corte_.shp")
shape
## [1] "corte_.shp"

Se guarda el shapefile con el que se va a realizar el corte en la variable “zona_corte”

zona_corte<-readOGR("C:\\Users\\Antonio\\Desktop\\PERCEPCION_REMOTA\\TALLER\\LABORATORIO 4\\shapefile_corte")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\Antonio\Desktop\PERCEPCION_REMOTA\TALLER\LABORATORIO 4\shapefile_corte", layer: "corte_"
## with 1 features
## It has 1 fields

Se verifica su extensión y se guarda en “e”

ext(zona_corte)
## SpatExtent : 415123.583251953, 550590.520996094, 4547057.90252686, 4649187.27337646 (xmin, xmax, ymin, ymax)
e<-ext(415123.583251953, 550590.520996094, 4547057.90252686, 4649187.27337646)
e
## SpatExtent : 415123.583251953, 550590.520996094, 4547057.90252686, 4649187.27337646 (xmin, xmax, ymin, ymax)

Se realiza el corte entre la imagen y las dimensiones del shapefile.

ndvi <- crop(ndvi, e)

se visualiza las caraceristricas del nuevo corte

ndvi
## class       : SpatRaster 
## dimensions  : 3404, 4516, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 415125, 550605, 4547055, 4649175  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs 
## source      : memory 
## name        :       NIR 
## min value   :        -1 
## max value   : 0.6598891

Se convierte el raster en un vector/ matriz

nr <- values(ndvi)

Se revisa su información:

str(ndvi)
## Formal class 'SpatRaster' [package "terra"] with 1 slot
##   ..@ ptr:Reference class 'Rcpp_SpatRaster' [package "terra"] with 16 fields
##   .. ..$ depth    : num 0
##   .. ..$ extent   :Reference class 'Rcpp_SpatExtent' [package "terra"] with 2 fields
##   .. .. ..$ valid : logi TRUE
##   .. .. ..$ vector: num [1:4] 415125 550605 4547055 4649175
##   .. .. ..and 27 methods, of which 13 are  possibly relevant:
##   .. .. ..  align, as.points, ceil, compare, finalize, floor, initialize,
##   .. .. ..  intersect, round, sample, sampleRandom, sampleRegular, union
##   .. ..$ filenames: chr ""
##   .. ..$ hasRange : logi TRUE
##   .. ..$ hasTime  : logi FALSE
##   .. ..$ hasValues: logi TRUE
##   .. ..$ inMemory : logi TRUE
##   .. ..$ messages :Reference class 'Rcpp_SpatMessages' [package "terra"] with 2 fields
##   .. .. ..$ has_error  : logi FALSE
##   .. .. ..$ has_warning: logi FALSE
##   .. .. ..and 18 methods, of which 4 are  possibly relevant:
##   .. .. ..  finalize, getError, getWarnings, initialize
##   .. ..$ names    : chr "NIR"
##   .. ..$ origin   : num [1:2] 15 15
##   .. ..$ range_max: num 0.66
##   .. ..$ range_min: num -1
##   .. ..$ res      : num [1:2] 30 30
##   .. ..$ time     : num 0
##   .. ..$ timestep : chr "seconds"
##   .. ..$ units    : chr ""
##   .. ..and 178 methods, of which 164 are  possibly relevant:
##   .. ..  addSource, adjacent, aggregate, align, apply, area_by_value,
##   .. ..  arith_numb, arith_rast, as_points, as_polygons, atan2, bilinearValues,
##   .. ..  boundaries, buffer, canProcessInMemory, cellFromRowCol,
##   .. ..  cellFromRowColCombine, cellFromXY, chunkSize, clamp, classify,
##   .. ..  colFromX, collapse_sources, combineSources, compare_geom,
##   .. ..  couldBeLonLat, count, cover, createAttributes, createCategories, crop,
##   .. ..  cum, deepcopy, dense_extent, disaggregate, expand, extCells,
##   .. ..  extractCell, extractVector, finalize, flip, focal, focalValues, freq,
##   .. ..  geometry, get_aggregate_dims, get_aggregates, get_crs,
##   .. ..  get_sourcenames, get_sourcenames_long, getAttributes, getBlockSize,
##   .. ..  getCategories, getColors, getError, getMessage, getNAflag,
##   .. ..  getRasterAtt, getValues, getWarnings, global, global_weighted_mean,
##   .. ..  gridDistance, has_error, has_warning, hasAttributes, hasCategories,
##   .. ..  hasColors, hasWindow, initf, initialize, initv, is_in, is_in_cells,
##   .. ..  isfinite, isGeographic, isGlobalLonLat, isinfinite, isLonLat, isnan,
##   .. ..  logic_numb, logic_rast, makeCategorical, mask_raster, mask_vector,
##   .. ..  math, math2, mem_needs, modal, ncol, nlyr, nlyrBySource, nrow, nsrc,
##   .. ..  polygonize, quantile, rapply, rappvals, rastDistance, rasterize,
##   .. ..  readAll, readStart, readStop, readValues, rectify, removeCategories,
##   .. ..  removeWindow, replace, reverse, rotate, rowColFromCell, rowFromY,
##   .. ..  rst_area, same, sampleRandomRaster, sampleRandomValues,
##   .. ..  sampleRegularRaster, sampleRegularValues, scale, selRange, separate,
##   .. ..  set_crs, set_depth, set_resolution, set_sourcenames,
##   .. ..  set_sourcenames_long, set_units, setAttributes, setCategories,
##   .. ..  setColors, setDepth, setNAflag, setNames, setRange, settime, setTime,
##   .. ..  setUnit, setValues, setWindow, shift, size, sources_to_disk, spatinit,
##   .. ..  stretch, subset, sum_area, summary, summary_numb, terrain, transpose,
##   .. ..  trig, trim, unique, vectCells, vectDistance, warp, writeRaster,
##   .. ..  writeStart, writeStop, writeValues, xFromCol, xyFromCell, yFromRow,
##   .. ..  zonal

Los values(ndvi) se convirtieron en ndviSpatRaster (matriz). Ahora se emplea kmeans que realiza un agrupamiento de los valores de la matriz

Se fija una semilla set.seed para la generación de números pseudoaleatorios

set.seed(99)

Se crean 10 clústeres que permiten 500 iteraciones, y comienzan con 5 conjuntos aleatorios usando el método "Lloyd

kmncluster <- kmeans(na.omit(nr), centers = 10, iter.max = 500, nstart = 5, algorithm="Lloyd")

Clase de kmncluster

str(kmncluster)
## List of 9
##  $ cluster     : int [1:15000034] 2 2 2 6 6 2 2 2 6 6 ...
##  $ centers     : num [1:10, 1] -0.0474 0.1626 0.5648 0.4581 0.2999 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:10] "1" "2" "3" "4" ...
##   .. ..$ : chr "NIR"
##  $ totss       : num 642827
##  $ withinss    : num [1:10] 242 386 430 408 397 ...
##  $ tot.withinss: num 3845
##  $ betweenss   : num 638982
##  $ size        : int [1:10] 3036646 708625 1006257 2191871 1294108 519930 1031298 1899749 1552269 1759281
##  $ iter        : int 333
##  $ ifault      : NULL
##  - attr(*, "class")= chr "kmeans"

Se utiliza el objeto ndvi para establecer los valores del clúster en un nuevo ráster

knr <- ndvi
values(knr) <- kmncluster$cluster
## Warning: [setValues] length of values does not match the number of cells
knr
## class       : SpatRaster 
## dimensions  : 3404, 4516, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 415125, 550605, 4547055, 4649175  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs 
## source      : memory 
## name        : NIR 
## min value   :   1 
## max value   :  10

El knr es un SpatRaster el cual creo gropos de 1 a 10 (maximo y mínimo), para poder determinar que cobertura corresponde se puede realizar mediante la ubicación de las imágenes (NDVI vs Kmeans) y definiendo un color para cada agrupación

mycolor <- c("#FF99FF","#0066FF","#00FF4D","#00FF66","#FF9900","#00ff00","#cbbeb5",
             "#FF0099","#CCFF00","#33FF00","#00FF66","#FF9900")
par(mfrow = c(1,2))
plot(ndvi, col = rev(terrain.colors(10)), main = "Landsat-NDVI")
plot(knr, main = 'Unsupervised classification', col=mycolor, type="classes")

Se realiza el corte a la imagen satelital.

landsat1 <- crop(landsat, e)
landsat1
## class       : SpatRaster 
## dimensions  : 3404, 4516, 6  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 415125, 550605, 4547055, 4649175  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=16 +datum=WGS84 +units=m +no_defs 
## source      : spat_LYK5BuCtv1oYuSd.tif 
## names       :  blue, green,   red,   NIR, SWIR1, SWIR2 
## min values  :     0,     0,     0,     0,     0,     0 
## max values  : 44760, 45410, 48261, 65535, 65535, 65535

Se plotean las tres imagenes

par(mfrow = c(2,2))
plot(knr, main = 'Unsupervised classification', col=mycolor, type="classes")
plot(ndvi, col = rev(terrain.colors(10)), main = "Landsat-NDVI")
plotRGB(landsat1, r=3, g=2, b=1, colNA="white", main="Landsat 8 Color Verdadero")

Ahora se etiqueta con el apoyo de la imagen en color verdadero para distiguir las coberturas visualmente.

eticq=c("Suelo descubierto","Edificaciones","Vegetación","Vegetación","Vegetación","Edificaciones","vegetación","Vegetación","Vegetación","Agua")

Este cuaderno se hizo basado en el del profesor Ivan Lizarazo que puede ser consultado en l siguiente link: https://rspatial.org/terra/rs/4-unsupclassification.html

sessionInfo()
## R version 4.0.4 (2021-02-15)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18363)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Spanish_Spain.1252  LC_CTYPE=Spanish_Spain.1252   
## [3] LC_MONETARY=Spanish_Spain.1252 LC_NUMERIC=C                  
## [5] LC_TIME=Spanish_Spain.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] corrplot_0.84 rgeos_0.5-5   terra_1.0-10  ggplot2_3.3.3 rgdal_1.5-23 
## [6] raster_3.4-5  sp_1.4-5     
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.6        highr_0.8         compiler_4.0.4    pillar_1.5.1     
##  [5] tools_4.0.4       digest_0.6.27     evaluate_0.14     lifecycle_1.0.0  
##  [9] tibble_3.1.0      gtable_0.3.0      lattice_0.20-41   pkgconfig_2.0.3  
## [13] rlang_0.4.10      DBI_1.1.1         yaml_2.2.1        xfun_0.21        
## [17] withr_2.4.1       stringr_1.4.0     dplyr_1.0.5       knitr_1.31       
## [21] generics_0.1.0    vctrs_0.3.6       grid_4.0.4        tidyselect_1.1.0 
## [25] glue_1.4.2        R6_2.5.0          fansi_0.4.2       rmarkdown_2.7    
## [29] purrr_0.3.4       magrittr_2.0.1    scales_1.1.1      codetools_0.2-18 
## [33] htmltools_0.5.1.1 ellipsis_0.3.1    assertthat_0.2.1  colorspace_2.0-0 
## [37] utf8_1.1.4        stringi_1.5.3     munsell_0.5.0     crayon_1.4.1