Realizar cálculos con datos raster

Muchas veces queremos combinar y realizar cálculos utilizando datos rasters, que como resultado tenemos un raster nuevo de salida.

La Función overlay crea un nuevo raster basado en dos o más rasters. Para utilizar overlay se debe crear una función para establecer la forma en que se van a combinar los rasters

Como ejemplo vamos a realizar el cálculo de un índice de vegetación (NDVI), utilizando las bandas 5 y 4 de una imagen landsat 8.

# Cargar librerias
library(raster)

## Loading required package: sp

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 3.4.4

#Establecer el directorio
setwd("D:\\maestria\\percepcion remota avanzada\\ejerciciosR")

#Subir bandas y asignamos una variable a cada una
ls8_b5 <- raster("LC08_L1TP_010062_20161120_20170318_01_T1_B5.tif")
ls8_b4 <- raster("LC08_L1TP_010062_20161120_20170318_01_T1_B4.tif")

#Visualizamos cada banda
par(mfrow=c(1,2))
plot(ls8_b5)
plot(ls8_b4)

Podemos calcular el ndvi de 2 formas, La primera realizando una operación matemática simple y la segunda utilizando la función overlay

#haciendo operaciones matematicas
ndvi_1 <- (ls8_b5-ls8_b4)/(ls8_b5+ls8_b4)

#utilizando la funcion overlay la cual crea un nuevo raster basado en 2 o mas rasters
ndvi_2 <- overlay(ls8_b5, ls8_b4, fun=function(x,y){(x-y)/(x+y)})

# Ploteamos ndvi_1 y ndvi_2 y observamos que los resultados fueron iguales 
par(mfrow=c(1,2))
plot(ndvi_1)
plot(ndvi_2)

Para conocer la distribución de los datos obtenidos de ndvi, generamos un histograma

#transorfar los datos de un RasterLayer a un dataframe
df_ndvi <- as.data.frame(ndvi_1, xy= TRUE)
df_ndvi <- na.omit(df_ndvi)

ggplot(df_ndvi) +
    geom_histogram(aes(layer), bins = 50)