En esta guía se comparte las pautas para realizar de preprocesamiento en las imágenes CBERS 4A, es decir la conversión de nivel digital a nivel de reflectancia al tope de la atmosfera (TOA), además, se muestran los procesos aritméticos para identificar objetos geográficos mediante análisis de índices espectrales, y luego los criterios de agrupamiento de pixeles en rangos de valores para asociar con los objetos geográficos de interés, el propósito de este paso a paso, es mostrar el potencial de Rstudio para el procesamiento de imágenes satelitales y así mismo, también mostrar el potencial de las imágenes CBERS 4A para evaluar y monitorear los recursos naturales que abarquen superficies extensas.
# cargando librerias
library(sp)
library(raster)
library(ggplot2)
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-1library(RStoolbox)
#E:/Imagenes CBERS4A/ParaR
setwd("E:/A_manuscrito/Para ponencia")
files <- list.files(pattern='.sdat')
files## [1] "Band 1.sdat" "Band 1.sdat.aux.xml" "Band 2.sdat"
## [4] "Band 2.sdat.aux.xml" "Band 3.sdat" "Band 3.sdat.aux.xml"
## [7] "Band 4.sdat" "Band 4.sdat.aux.xml" "PAN.sdat"
## [10] "PAN.sdat.aux.xml"b1 <- raster("E:/A_manuscrito/Para ponencia/Band 1.sdat")
b2 <- raster("E:/A_manuscrito/Para ponencia/Band 2.sdat")
b3 <- raster("E:/A_manuscrito/Para ponencia/Band 3.sdat")
b4 <- raster("E:/A_manuscrito/Para ponencia/Band 4.sdat")
#Características de los datos de las bandas multiespectrales, tamaño de pixel 8 metros.
b1## class : RasterLayer
## dimensions : 1355, 1642, 2224910 (nrow, ncol, ncell)
## resolution : 8, 8 (x, y)
## extent : 338258, 351394, 8745258, 8756098 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : Band 1.sdat
## names : Band_1
## values : -32768, 32767 (min, max)# Conversion de nivel digital a radiancia. ver articulo https://doi.org/10.1016/j.rse.2018.09.017
rad <- function(x1,x2,x3,x4){
rb1 <- 1.68*x1+(0)
rb2 <- 1.62*x2+(0)
rb3 <- 1.59*x3+(0)
rb4 <- 1.24*x4+(0)
return(list(rb1,rb2,rb3,rb4))
}
rabs<- stack(rad(b1,b2,b3,b4))
#Cálculo de distancia sol a tierra unidades astronómicas.
dist<- function(j=228){
da <- (1-0.01672*cos(0.9856*(j-4)*(pi/180)))
return(da)}# ingresamos el valor del dia juliano en este caso es 228
d<-dist()
d## [1] 1.012662#Cálculo de angulo cenital.
ac <- function(t=58.4011){
a <- cos((pi/180)*(90-t))
return(a)
}
z<- ac()
z## [1] 0.851737# calculo de reflectancia
#para determinar valores de Esun https://doi.org/10.1016/j.rse.2018.09.017
ref <- function(){
rf1<- ((pi)*rabs$Band_1*d^2)/(1958*z)
rf2<- ((pi)*rabs$Band_2*d^2)/(1852*z)
rf3<- ((pi)*rabs$Band_3*d^2)/(1559*z)
rf4<- ((pi)*rabs$Band_4*d^2)/(1091*z)
return(list(rf1,rf2,rf3,rf4))
}
rfbs<-stack(ref())
rfbs## class : RasterStack
## dimensions : 1355, 1642, 2224910, 4 (nrow, ncol, ncell, nlayers)
## resolution : 8, 8 (x, y)
## extent : 338258, 351394, 8745258, 8756098 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : Band_1, Band_2, Band_3, Band_4
## min values : 0.2791054, 0.2051346, 0.2044561, 0.2106523
## max values : 3.320057, 3.384721, 3.946389, 4.397905# generando el proceso de rescalamineto
res <- function(){
res1<-((rfbs$Band_1-0.2791054)/(3.320057-0.2791054))
res2<-((rfbs$Band_2-0.2051346)/(3.384721-0.2051346))
res3<-((rfbs$Band_3-0.2044561)/(3.946389-0.2044561))
res4<-((rfbs$Band_4-0.2106523)/(4.397905-0.2106523))
return(list(res1,res2,res3,res4))
}
resc<-stack(res())
resc## class : RasterStack
## dimensions : 1355, 1642, 2224910, 4 (nrow, ncol, ncell, nlayers)
## resolution : 8, 8 (x, y)
## extent : 338258, 351394, 8745258, 8756098 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : Band_1, Band_2, Band_3, Band_4
## min values : 1.185213e-08, 7.796832e-09, 7.050049e-09, 6.829703e-09
## max values : 0.9999998, 1.0000001, 1.0000000, 0.9999999# Explorando los datos preprocesados
library(ggplot2)
library(RStoolbox)
summary(resc$Band_1)## Band_1
## Min. 1.185213e-08
## 1st Qu. 1.056563e-01
## Median 1.259338e-01
## 3rd Qu. 1.579509e-01
## Max. 9.999998e-01
## NA's 1.961200e+04hist(resc$Band_4)
plot(resc$Band_2)
#Combinacion de bandas.
ggRGB(resc,r=3,g=2,b=1,stretch = "lin") + ggtitle("Combinacion bandas 321") 
# Haciendo los calculos de inidiceS de vegetacion.
#Calculo de NDVI
index <- function(x,y){
ndvi <- (x - y ) / (y + x)
return(ndvi)
}
NDVI<-index(resc$Band_4,resc$Band_3)
hist(NDVI)
plot(NDVI, col = rev(terrain.colors(30)), main = 'NDVI con imagenes CBERS 4A')
#Calculo de NDWI
index <- function(x,y){
ndvi <- (x - y ) / (y + x)
return(ndvi)
}
NDWI <-index(resc$Band_2,resc$Band_4)
hist(NDWI)
plot(NDWI, col = rev(terrain.colors(30)), main = 'NDWI con imagenes CBERS 4A')
#Calculo de EVI
index <- function(x,y,Z){
ndvi <- (x - y ) / (y +6*(x)-7.5*(z)+1)
return(ndvi)
}
EVI <-index(rfbs$Band_4,rfbs$Band_3,rfbs$Band_1)
hist(EVI)
#Agrupando pixeles.
veg1 <- calc(EVI, function(x){x[x < -0.3] <- NA; return(x)})
veg1 <- calc(veg1, function(x){x[x > 0.4] <- NA; return(x)})
plot(veg1, main = 'EVI con imagenes CBERS 4A', col = rev(terrain.colors(12)))
#Calculo de SAVI
index <- function(x,y){
ndvi <- (x - y ) / ((y + x)*(1.5))
return(ndvi)
}
SAVI <-index(resc$Band_4,resc$Band_3)
hist(SAVI)
plot(SAVI, col = rev(terrain.colors(30)), main = 'SAVI con imagenes CBERS 4A')
#Calculo de TCG,ese ejecuto en refencia al articulo doi:10.1631/jzus.B1100088
index <- function(w,x,y,z){
ndvi <- (-0.494*(w)-0.318*(x)-0.324*(y)+0.741*(z))
return(ndvi)
}
TCG <-index(resc$Band_1,resc$Band_2,resc$Band_3,resc$Band_4)
hist(TCG)
plot(TCG, col = rev(terrain.colors(30)), main = 'SAVI con imagenes CBERS 4A')
Despues de encontrar las nuevas bandas del produnto de las operaciones aritmeticas se salvaron estas nuevas bandas en el archivo de disco local, y para continuar nuevamente se importaron para continuar con el proceso de agrupamiento de pixeles.
vegr <- reclassify(TCG, c(-Inf,-0.18,NA, -0.18,0,NA, 0,0.16,NA, 0.16,Inf,4))
plot(vegr,col = rev(terrain.colors(5)), main = 'Clasificación coberturas TCG')
#
vegr1 <- reclassify(TCG,c(-Inf,-0.18,1, -0.18,0,2, 0,0.16,3, 0.16,Inf,4))
plot(vegr1,col = rev(terrain.colors(5)), main = 'Clasificación coberturas TCG')
# converientiendo valores de raster en vector
r.to.poly<-rasterToPolygons(vegr, dissolve = T)
# definiendo su projecccion
proj4string(TCG) <- proj4string(r.to.poly)
# mostrando dato vectorial
plot(r.to.poly,axe=TRUE)