5. Trabajar con raster multi-bandas en R

5 Trabajar con raster multi-bandas en R

Este codigo fue tomado de los tutoriales publicados por Data Carpentry, realizado por Leah A. Wasser, Megan A. Jones, Zack Brym, Kristina Riemer, Jason Williams, Jeff Hollister, Mike Smorul, Joseph Stachelek como parate del curso Introduccion a los Datos Geoespaciales con R

Introduccion

Este tutorial explora como importar y plotear raster multi-banda en R.

• Cargar librerias

library(raster)

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

## Loading required package: sp

library(rgdal)

## rgdal: version: 1.2-16, (SVN revision 701)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.0, released 2017/04/28
##  Path to GDAL shared files: C:/Users/Natalia/Documents/R/win-library/3.4/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/Natalia/Documents/R/win-library/3.4/rgdal/proj
##  Linking to sp version: 1.2-7

library(sp)
library(ggplot2)

Para leer un objeto Raster Stack en R se utiliza la funcion raster(), esta lee unicamente la primera banda

Stack_zona <-
  raster("C:\\Maestria en Geomatica\\Percepcion remota avanzada\\Talleres\\Stack_zonaestudio.tif")

setwd("C:\\Maestria en Geomatica\\Percepcion remota avanzada\\Talleres")
getwd()

## [1] "C:/Maestria en Geomatica/Percepcion remota avanzada/Talleres"

GDALinfo("Stack_zonaestudio.tif")

## Warning in GDALinfo("Stack_zonaestudio.tif"): statistics not supported by
## this driver

## Warning in GDALinfo("Stack_zonaestudio.tif"): statistics not supported by
## this driver

## Warning in GDALinfo("Stack_zonaestudio.tif"): statistics not supported by
## this driver

## rows        2993 
## columns     3283 
## bands       6 
## lower left origin.x        625395 
## lower left origin.y        397485 
## res.x       30 
## res.y       30 
## ysign       -1 
## oblique.x   0 
## oblique.y   0 
## driver      GTiff 
## projection  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## file        Stack_zonaestudio.tif 
## apparent band summary:
##    GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float32          FALSE           0          1       3283
## 2 Float32          FALSE           0          1       3283
## 3 Float32          FALSE           0          1       3283
## 4 Float32          FALSE           0          1       3283
## 5 Float32          FALSE           0          1       3283
## 6 Float32          FALSE           0          1       3283
## apparent band statistics:
##          Bmin       Bmax    Bmean      Bsd
## 1        7366      38151 9161.524 1351.857
## 2        6379      38463 8876.713 1416.138
## 3        5780      40551 7997.838 1740.965
## 4 -4294967295 4294967295       NA       NA
## 5 -4294967295 4294967295       NA       NA
## 6 -4294967295 4294967295       NA       NA
## Metadata:
## AREA_OR_POINT=Area

Se necesita convertir este dato en data frame para podrer realizar el ploteo con ggplot

Stack_zona_df  <- raster::as.data.frame(Stack_zona, xy = TRUE)
Stack_zona_df <- na.omit(Stack_zona_df)

ggplot() +
 geom_raster(data = Stack_zona_df , aes(x = x, y = y, alpha = Stack_zonaestudio))

Valores de los datos de imagen raster

Los raster contienen valores entre 10000 y 40000 . Esos valores representan grados brightness asociados con la banda de la imagen. En el caso de una imagen RGB (rojo, verde y azul), la banda 1 es la roja. Cuando se plotea la banda roja, los numeros grandes representan pixeles con mas rojo en ellos (un fuerte reflejo del rojo). Numeros pequeños representan pixeles con menos rojo (bajo reflejo del rojo). Para plotear una imagen RGB, es necesario mezclar los valores de rojo + verde + azul en una solo color para crear una imagen a color completa - un color de imagen similar al que crean las camaras digitales.

Importar una banda especifica

Par importar una banda especifica se utiliza la funcion raster(), especificando que banda se quiere band = N (N representa el numero de la banda con la cual se desea trabajar). Ej. si se desea importar la banda verde y esta es la numero 2 se utiliza band = 2

Band2 <-
  raster("C:\\Maestria en Geomatica\\Percepcion remota avanzada\\Talleres\\Stack_zonaestudio.tif",
           band = 2)

Estos datos pueden convertirse a data frame y plotear la banda seleccionada

Band2_df  <- raster::as.data.frame(Band2, xy = TRUE)
Band2_df <- na.omit(Band2_df)

ggplot() +
 geom_raster(data = Band2_df , aes(x = x, y = y, alpha = Stack_zonaestudio))

Raster Stacks en R

Para trabajar con todas las bandas de la imagen se utiliza como objeto el RasterStack. Despues se pueden plotear diferentes composiciones de la imagen.

Para traer todas las bandas como un raster multi-banda, se utiliza la funcion stack().

stack_zona <-
  stack("C:\\Maestria en Geomatica\\Percepcion remota avanzada\\Talleres\\Stack_zonaestudio.tif")

Ver los atributos del objeto stack

stack_zona

## class       : RasterStack 
## dimensions  : 2993, 3283, 9826019, 6  (nrow, ncol, ncell, nlayers)
## resolution  : 30, 30  (x, y)
## extent      : 625395, 723885, 397485, 487275  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names       : Stack_zonaestudio.1, Stack_zonaestudio.2, Stack_zonaestudio.3, Stack_zonaestudio.4, Stack_zonaestudio.5, Stack_zonaestudio.6 
## min values  :                7366,                6379,                5780,                   ?,                   ?,                   ? 
## max values  :               38151,               38463,               40551,                   ?,                   ?,                   ?

Ver los atributos de cada banda de forma independiente

stack_zona@layers

## [[1]]
## class       : RasterLayer 
## band        : 1  (of  6  bands)
## dimensions  : 2993, 3283, 9826019  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 625395, 723885, 397485, 487275  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : C:\Maestria en Geomatica\Percepcion remota avanzada\Talleres\Stack_zonaestudio.tif 
## names       : Stack_zonaestudio.1 
## values      : 7366, 38151  (min, max)
## 
## 
## [[2]]
## class       : RasterLayer 
## band        : 2  (of  6  bands)
## dimensions  : 2993, 3283, 9826019  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 625395, 723885, 397485, 487275  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : C:\Maestria en Geomatica\Percepcion remota avanzada\Talleres\Stack_zonaestudio.tif 
## names       : Stack_zonaestudio.2 
## values      : 6379, 38463  (min, max)
## 
## 
## [[3]]
## class       : RasterLayer 
## band        : 3  (of  6  bands)
## dimensions  : 2993, 3283, 9826019  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 625395, 723885, 397485, 487275  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : C:\Maestria en Geomatica\Percepcion remota avanzada\Talleres\Stack_zonaestudio.tif 
## names       : Stack_zonaestudio.3 
## values      : 5780, 40551  (min, max)
## 
## 
## [[4]]
## class       : RasterLayer 
## band        : 4  (of  6  bands)
## dimensions  : 2993, 3283, 9826019  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 625395, 723885, 397485, 487275  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : C:\Maestria en Geomatica\Percepcion remota avanzada\Talleres\Stack_zonaestudio.tif 
## names       : Stack_zonaestudio.4 
## values      : NA, NA  (min, max)
## 
## 
## [[5]]
## class       : RasterLayer 
## band        : 5  (of  6  bands)
## dimensions  : 2993, 3283, 9826019  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 625395, 723885, 397485, 487275  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : C:\Maestria en Geomatica\Percepcion remota avanzada\Talleres\Stack_zonaestudio.tif 
## names       : Stack_zonaestudio.5 
## values      : NA, NA  (min, max)
## 
## 
## [[6]]
## class       : RasterLayer 
## band        : 6  (of  6  bands)
## dimensions  : 2993, 3283, 9826019  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 625395, 723885, 397485, 487275  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : C:\Maestria en Geomatica\Percepcion remota avanzada\Talleres\Stack_zonaestudio.tif 
## names       : Stack_zonaestudio.6 
## values      : NA, NA  (min, max)

Si se tienen muchas bandas, es posible especificar las bandas de las cuales quiero ver los atributos utilizando un indice de valor

stack_zona[[2]]

## class       : RasterLayer 
## band        : 2  (of  6  bands)
## dimensions  : 2993, 3283, 9826019  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 625395, 723885, 397485, 487275  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : C:\Maestria en Geomatica\Percepcion remota avanzada\Talleres\Stack_zonaestudio.tif 
## names       : Stack_zonaestudio.2 
## values      : 6379, 38463  (min, max)

Es posible tambien utilizar la función ggplot para plotear cualquier layer o el RasterStack. Es importante recordar que primero es necesario convertirlo a un data frame

stack_zona_df <- raster::as.data.frame(stack_zona, xy = TRUE)
stack_zona_df <- na.omit(stack_zona_df)

Cada banda en el RasterStack toma su propia columna en el data frame.

str(stack_zona_df)

## 'data.frame':    5110153 obs. of  8 variables:
##  $ x                  : num  653820 653790 653820 653790 653820 ...
##  $ y                  : num  487260 487230 487230 487200 487200 ...
##  $ Stack_zonaestudio.1: num  9452 10517 9551 8947 8963 ...
##  $ Stack_zonaestudio.2: num  9055 9796 8920 8318 8432 ...
##  $ Stack_zonaestudio.3: num  8280 9203 8229 7579 7590 ...
##  $ Stack_zonaestudio.4: num  18623 18108 17587 17353 17709 ...
##  $ Stack_zonaestudio.5: num  11613 12215 11265 10654 10820 ...
##  $ Stack_zonaestudio.6: num  8684 9746 8603 7949 7933 ...
##  - attr(*, "na.action")=Class 'omit'  Named int [1:4715866] 1 2 3 4 5 6 7 8 9 10 ...
##   .. ..- attr(*, "names")= chr [1:4715866] "1" "2" "3" "4" ...

• Histograma de la primera banda

ggplot() +
  geom_histogram(data = stack_zona_df, aes(Stack_zonaestudio.1))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

• Histograma con ajustes para la visualización

ggplot() +
  geom_histogram(data = stack_zona_df, aes(Stack_zonaestudio.1)) +
   scale_x_continuous(limits = c(8000,14000))

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## Warning: Removed 72179 rows containing non-finite values (stat_bin).

## Warning: Removed 1 rows containing missing values (geom_bar).

Ploteo de la segunda banda del raster

ggplot() +
  geom_raster(data = stack_zona_df, aes(x = x, y = y, alpha = Stack_zonaestudio.2))

Crear una Imagen con 3 bandas

Para hacer una imagen con tres a color en R se utiliza la funcion plotRGB()

Esta funcion permite:

1. Identificar que bandas se quiere utilizar en las regiones roja, verde y azul. La funcion plotRGB() predetermina el orden de las bandas asi: 1 = rojo, 2 = verde y 3 = azul. Sinembargo, se pueden definir las bandas que se desean plotear manualmente. La definicion manual de las bandas es util si se tiene por ejemplo la banda infrarojo cercano y se desea crear una imagen con color infrarojo.
2. Ajustar el strech de una imagen permite incrementar o disminuir el contraste.

Como se mencionó anteriormente para plotear la imagen de 3 bandas se utiliza la funcion plotRGB() directamente sobre el objeto RasterStack (no se necesita el dataframe ya que esta funcion no es parte del paquete ggplot)

Modificaciones de la imagen cuando se utiliza el stretch, puede mejorar la claridad y contraste de la imagen utilizando stretch=“lin” o stretch=“hist”

Cuando el rango de los valores de brillo del pixel estan cera a 0, la imagen es oscura por defecto. Es posible realizar stretch (estiramiento) de valores en un rango completo en el rango de los valores, para aumentar el contraste visual de la imagen.

plotRGB(stack_zona,
        r = 3, g = 2, b = 1,
        scale = 800,
        stretch = "lin")

plotRGB(stack_zona,
        r = 3, g = 2, b = 1,
        scale = 800,
        stretch = "hist")

En este caso el estiramiento cambia notoriamente el contraste de la imagen.

Raster Stack vs Raster Brick en R

Los objetos RasterStack y RasterBrik ambos almacenan multiples bandas. Sinembargo, la forma como almacenan las bandas es diferente. Las bandas en RasterStack se almacenan como enlaces al dato raster que se encuentra en algun lugar del computador. Un RasterBrick contiene todos los objetos guardados en el objeto actual en R. En muchos casos, se puede trabajar con RAsterBrick en la misma forma como se trabaja con RasterStack. Sinembargo un RasterBrick es generalmente mas eficiente y rapido de procesar - esto es importanto cuando se trabajan archivos muy pesados.

Es posible pasar de RasterStack a RasterBrick en R usando brick(nombre del stack). Se puede utilizar la funcion object.size() para comparar el RasterStack y el RasterBrick.

• Revision de tamaño del RasterStack

object.size(stack_zona)

## 78304 bytes

• Crear RasterBrick

brick_zona <- brick(stack_zona)

object.size(brick_zona)

## 471662400 bytes

En el RasterBrick todas las bandas estan almacenadas en el objeto actual. Por esto, el tamaño del RasterBrick es mucho mas grande que la del RasterStack

Para plotear los RasterBrick se utiliza tambien la funcion plotRGB()

plotRGB(brick_zona)