Preprocesamiento de imagenes satelitales y aplicación de indice de vegetación

Programa de ingeniería en topografía - Unidades Tecnologicas de Santander

Este proceso intenta mostrar las pautas de preprocesamiento de las imágenes digitales, la conversión de Nivel Digital (ND) a valores de radiancia, de valores a radiancia a valores de reflectancia, análisis de índice de vegetación utilizando algunos algoritmos, en R Markdonw en línea, siendo este proceso de suma importancia para el procesamiento y análisis espectral de los objetos mediante las imágenes digitales. Para realizar este proceso se utilizaron imágenes Landsat 8, que se encuentra disponible para el acceso en la página de Earth Explorer, así también, para la transformación y aplicación de las imágenes digitales han sido desarrolladas en el programa R Package y RStudio utilizando librerías básicas para el uso de los datos espaciales tales como: Raster, Rgdal, Sp, Rtoolbox, y sobre todo se intenta aplicar funciones para reducir el error en procesamiento de los datos, y de esta manera el procesamiento sea más amigable, considerando estos criterios, el propósito de esta guía es para mostrar el proceso y el uso de los códigos para el prerocesamiento de imágenes digitales y análisis de índices de diferencia normalizada de vegetación, a los aprendices de los programas de ingeniería y de las ciencias de la tierra y sociedad civil interesada, teniendo en cuenta el ámbito regional nacional e internacional.

library(sp)
library(raster)
library(rgeos)

## rgeos version: 0.5-2, (SVN revision 621)
##  GEOS runtime version: 3.6.1-CAPI-1.10.1 
##  Linking to sp version: 1.3-2 
##  Polygon checking: TRUE

library(rgdal)

## rgdal: version: 1.4-8, (SVN revision 845)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/Usuario/Documents/R/win-library/3.6/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/Usuario/Documents/R/win-library/3.6/rgdal/proj
##  Linking to sp version: 1.4-1

# importamos los datos raster 
#E:\Datos vilca\sancarlos
setwd("E:/Datos vilca/sancarlos")
files <- list.files(pattern='.tif')
files

##  [1] "B1.tif"           "B1.tif.enp"       "B2.tif"           "B2.tif.enp"      
##  [5] "B3.tif"           "B3.tif.enp"       "B4.tif"           "B4.tif.enp"      
##  [9] "B5.tif"           "B5.tif.enp"       "B6.tif"           "B6.tif.enp"      
## [13] "B7.tif"           "B7.tif.enp"       "B8.tif"           "B8.tif.enp"      
## [17] "data.tif"         "data.tif.aux.xml" "data.tif.enp"     "data.tif.ovr"

b1 <- raster("E:/Datos vilca/sancarlos/B1.tif")
b2 <- raster("E:/Datos vilca/sancarlos/B2.tif")
b3 <- raster("E:/Datos vilca/sancarlos/B3.tif")
b4 <- raster("E:/Datos vilca/sancarlos/B4.tif")
b5 <- raster("E:/Datos vilca/sancarlos/B5.tif")
b6 <- raster("E:/Datos vilca/sancarlos/B6.tif")
b7 <- raster("E:/Datos vilca/sancarlos/B7.tif")

b4

## class      : RasterLayer 
## dimensions : 231, 223, 51513  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 475365, 482055, -1387305, -1380375  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : E:/Datos vilca/sancarlos/B4.tif 
## names      : B4 
## values     : -2147483648, 2147483647  (min, max)

plot(b4)

# Conversion de nivel digital a radiancia. 
rad <- function(x1,x2,x3,x4,x5,x6,x7){ 
  rb1 <- 0.012199*x1+ (-60.99626)
  rb2 <- 0.012492*x2+(-62.46091)
  rb3 <- 0.011511*x3+(-57.55720)
  rb4 <- 0.009707*x4+(-48.53547)
  rb5 <- 0.00594*x5+(-29.70129)
  rb6 <- 0.001477*x6+(-7.38644)
  rb7 <- 0.0004979*x7+(-2.48963)
  return(list(rb1,rb2,rb3,rb4,rb5,rb6,rb7))
}
rabs<- stack(rad(b1,b2,b3,b4,b5,b6,b7))

#Cálculo de distancia sol a tierra unidades astronómicas.
dist<- function(j=217){
  da <- (1-0.01672*cos(0.9856*(j-4)*(pi/180)))
  return(da)}# ingresamos el valor del dia juliano en este caso es 217
d<-dist()
d

## [1] 1.01449

#Cálculo de angulo cenital.
ac <- function(t=47.612){
  a <- cos((pi/180)*(90-t))
  return(a)
}
z<- ac()
z

## [1] 0.7385965

# calculo de reflectancia 
#para determinar valores de Esun http://www.gisagmaps.com/landsat-8-atco/
ref <- function(){
  rf2<- ((pi)*rabs$B2*d^2)/(2067*z)
  rf3<- ((pi)*rabs$B3*d^2)/(1893*z)
  rf4<- ((pi)*rabs$B4*d^2)/(1603*z)
  rf5<- ((pi)*rabs$B5*d^2)/(972.6*z)
  rf6<- ((pi)*rabs$B6*d^2)/(245*z)
  rf7<- ((pi)*rabs$B7*d^2)/(79.72*z)
  return(list(rf2,rf3,rf4,rf5,rf6,rf7))
}

rfbs<-stack(ref())
rfbs

## class      : RasterStack 
## dimensions : 231, 223, 51513, 6  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 475365, 482055, -1387305, -1380375  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names      :          B2,          B3,          B4,          B5,          B6,          B7 
## min values : 0.052751972, 0.030926804, 0.020171883, 0.014992871, 0.009422172, 0.006609358 
## max values :   0.2278663,   0.2681068,   0.3182099,   0.5108582,   0.6290514,   0.5327576

library(ggplot2)
library(RStoolbox)
summary(rfbs$B2)

##                 B2
## Min.    0.05275197
## 1st Qu. 0.08613985
## Median  0.09923572
## 3rd Qu. 0.11061194
## Max.    0.22786635
## NA's    0.00000000

#mostrando los datos raster en el mapa.
plot(rfbs$B2)

#Combinacion de bandas.
ggRGB(rfbs,r=3,g=2,b=1,stretch = "lin") + ggtitle("Combinacion bandas 321") 

# Haciendo los calculos de inidice de vegetacion.
index <- function(x,y){
  ndvi <- (y - x ) / (y + x) 
  return(ndvi)
}
NDVI<-index(rfbs$B3,rfbs$B4)

hist(NDVI)

plot(NDVI, col = rev(terrain.colors(30)), main = 'NDVI con imagenes Landsat 8')

# reclasificando en base a umbrales 
vegr <- reclassify(NDVI, c(-0.2,-0.1,1, -0.1,0,2, 0,0.05,3, 0.05,0.1,4, 0.1,0.15,5,0.15,Inf,6))
plot(vegr,col = rev(terrain.colors(5)), main = 'NDVI base a umbrales')

vegr1 <- reclassify(NDVI, c(-0.2,-0.1,NA, -0.1,0,NA, 0,0.05,NA, 0.05,0.1,NA, 0.1,0.15,5,0.15,Inf,6))
plot(vegr1,col = rev(terrain.colors(5)), main = 'NDVI base a umbrales')

#Agrupando pixeles.
veg1 <- calc(NDVI, function(x){x[x < -0.3] <- NA; return(x)})
veg1 <- calc(veg1, function(x){x[x > 0.3] <- NA; return(x)})
plot(veg1, main = 'Coberturas', col = rev(terrain.colors(8)))

#Agrupando Pixeles en rango de valores
veg2 <- calc(NDVI, function(x){x[x < 0.12] <- NA; return(x)})
veg2 <- calc(veg2, function(x){x[x > 0.4] <- NA; return(x)})
plot(veg2, main = 'Indicios de vegetacion \nNDVI > 0.12', col = rev(topo.colors(4)))

# superponiendo capas 

plotRGB(rfbs, r=4, g=3, b=2, axes=TRUE, stretch="lin", main="Landsat")
plot(veg2, add=TRUE, legend=FALSE)

#
plotRGB(rfbs, r=4, g=3, b=2, axes=FALSE, stretch="lin", main="Landsat")
plot(vegr1, add=TRUE, legend=FALSE)