librerias

library(terra)
library(raster)
library(sp)

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) 
writeRaster(landsat, filename="math.tif", overwrite=TRUE)
LSbrick <- brick(input('math',type = '.tif'))

Funcion Vegetation Index

vi <- function(img, k, i) {
  bk <- img[[k]]
  bi <- img[[i]]
  vi <- (bk - bi) / (bk + bi)
  return(vi)
}
ndvi <- vi(landsat,5,4)
plot(ndvi,col=rev(terrain.colors(100)),main='landsat NDVI')

Pruebas Indices de agua

reference:

NDWIG<- vi(LSbrick,3,5)
NDWIR<- vi(LSbrick,4,6)
MNDWI1<- vi(LSbrick,3,6)
MNDWI2<- vi(LSbrick,3,7)
AWEIsh<-(LSbrick$blue+(2.5*LSbrick$green)-(1.5*(LSbrick$NIR+LSbrick$SWIR1)) -(0.25*LSbrick$SWIR2))
AWEInsh<- (4*(LSbrick$green-LSbrick$SWIR1))-((0.25*LSbrick$NIR)+(2.75*LSbrick$SWIR1))

##Normalizar indices

norma <- function(img) {
  mx <- maxValue(img)
  mn <- minValue(img)
  norma <- 100-((100*(mx-img))/(mx-mn))
  return(norma)
}
nNDWIG=norma(NDWIG)
nMNDWI1=norma(MNDWI1)
nNDWIR=norma(NDWIR)
nMNDWI2=norma(MNDWI2)
nAWEIsh=norma(AWEIsh)
nAWEInsh=norma(AWEInsh)
plot(nNDWIG,col=rev(terrain.colors(100)),axes=F,main='landsat NDWIG')

plot(nNDWIR,col=rev(terrain.colors(100)),axes=F,main='landsat NDWIR')

plot(nMNDWI1,col=rev(terrain.colors(100)),axes=F,main='landsat MNDWI1')

plot(nMNDWI2,col=rev(terrain.colors(100)),axes=F,main='landsat MNDWI2')

plot(nAWEIsh,col=rev(terrain.colors(100)),axes=F,main='landsat AWEIsh')

plot(nAWEInsh,col=rev(terrain.colors(100)),axes=F,main='landsat AWEInsh')

histogramas

hist(ndvi,main = "Distribucion NDVI",xlab = "NDVI",ylab= "Frecuencia",col = rev(terrain.colors(40)),xlim = c(-0.25, 0.6), breaks = 30, xaxt = 'n')
axis(side=1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))

par(mfrow=c(2,2))
hist(nNDWIG,main="Distribucion NDWIg",xlab="NDVWIg",ylab="Frecuencia",col = rev(terrain.colors(45)), xlim = c(0, 100),ylim=c(0,20000), breaks = 40, xaxt = 'n')
axis(side=1, at = seq(0,100,10), labels = seq(0,100,10))
hist(nMNDWI2,main="Distribucion MNDWI2",xlab="MNDWI2",ylab="Frecuencia",col = rev(terrain.colors(45)), xlim = c(0, 100),ylim=c(0,20000), breaks = 40, xaxt = 'n')
axis(side=1, at = seq(0,100,10), labels = seq(0,100,10))
plot(nNDWIG,col=rev(terrain.colors(100)),axes=F,main='landsat NDWIG')
plot(nMNDWI2,col=rev(terrain.colors(100)),axes=F,main='landsat MNDWI2')

Umbrales, Thresholding

veg <- classify(ndvi, cbind(-Inf, 0.4, NA))
plot(veg, main='Vegetation')

Categorias

ranges <- c(-Inf,0.25,1,0.25,0.3,2,0.3,0.4,3,0.4,0.5,4,0.5,Inf, 5)
Mranges <- matrix(ranges, ncol = 3, byrow = T)
clasifi <- classify(ndvi, Mranges)
plot(clasifi,col = c('red','orange','yellow','green','darkgreen'),
     main = 'Umbrales de NDVI',legend=F)
legend(x ='topright',legend=c('1','2','3','4','5'), fill=c('red','orange','yellow','green','darkgreen'))

deteccion de cuerpos de agua con NDVI

rangeW <- c(-Inf,0.05,1,-1,Inf,2)
MrangeW <- matrix(rangeW, ncol=3, byrow=T)
clasifiW <- classify(ndvi,MrangeW)
plot(clasifiW,col = c('black','white'),main = 'Cuerpos de agua con NDVI',legend=F)
legend(x ='topright',legend=c('Agua'), fill=c('black'))

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] raster_3.3-7 sp_1.4-2     terra_1.0-10
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.5        codetools_0.2-18  lattice_0.20-41   digest_0.6.25    
##  [5] grid_4.0.4        magrittr_1.5      evaluate_0.14     highr_0.8        
##  [9] rlang_0.4.10      stringi_1.4.6     rmarkdown_2.7     rgdal_1.5-23     
## [13] tools_4.0.4       stringr_1.4.0     xfun_0.21         yaml_2.2.1       
## [17] compiler_4.0.4    htmltools_0.5.1.1 knitr_1.31