4. Calculos de Rasters 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

A menudo se necesitan combinar valores y realizar calculos entre rasters para crear nuevos rasters de salida. Esto incluye como restar un raster de otro usando matematica basica entre rasters y la funcion overlay(). Tambien incluye como extraerlos valores de pixeles de un conjunto de ubicaciones, por ejemplo el buffer de una region alrededor alrededor de la ubicacion delas parcelas de un sitio en campo.

  • 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)
  • Importar datos Raster DEM y DEM hillshade
DEM <- raster("C:\\Maestria en Geomatica\\Percepcion remota avanzada\\Talleres\\DEM_zona.tif")

# import DTM hillshade
Hillshade <- raster("C:\\Maestria en Geomatica\\Percepcion remota avanzada\\Talleres\\Hillshade_DEM1.tif")

# Imagen B2
Imagen <- raster("C:\\Maestria en Geomatica\\Percepcion remota avanzada\\Talleres\\Landsat_recorte_B2.tif")
  • Convertir en data frame
DEM_zona_df <- raster::as.data.frame(DEM, xy = TRUE)
DEM_zona_df <- na.omit(DEM_zona_df)
str(DEM_zona_df)
## 'data.frame':    553959 obs. of  3 variables:
##  $ x       : num  1051460 1051460 1051369 1051460 1051369 ...
##  $ y       : num  979087 978996 978905 978905 978814 ...
##  $ DEM_zona: int  2956 2839 2706 2690 2638 2633 2604 2560 2560 2540 ...
##  - attr(*, "na.action")=Class 'omit'  Named int [1:514069] 1 2 3 4 5 6 7 8 9 10 ...
##   .. ..- attr(*, "names")= chr [1:514069] "1" "2" "3" "4" ...
Hill_df <- raster::as.data.frame(Hillshade, xy = TRUE)
Hill_df <- na.omit(Hill_df)
str(Hill_df)
## 'data.frame':    549211 obs. of  3 variables:
##  $ x             : num  1051460 1051369 1051460 1051278 1051369 ...
##  $ y             : num  978723 978632 978632 978541 978541 ...
##  $ Hillshade_DEM1: int  7 98 48 111 112 56 67 63 12 6 ...
##  - attr(*, "na.action")=Class 'omit'  Named int [1:518817] 1 2 3 4 5 6 7 8 9 10 ...
##   .. ..- attr(*, "names")= chr [1:518817] "1" "2" "3" "4" ...
Imagen_df <- raster::as.data.frame(Imagen, xy = TRUE)
Imagen_df <- na.omit(Imagen_df)
str(Imagen_df)
## 'data.frame':    5012512 obs. of  3 variables:
##  $ x                 : num  635820 635850 635880 635910 635940 ...
##  $ y                 : num  476040 476040 476040 476040 476040 ...
##  $ Landsat_recorte_B2: num  7796 7920 7951 8034 8240 ...
##  - attr(*, "na.action")=Class 'omit'  Named int [1:4336366] 1 2 3 4 5 6 7 8 9 10 ...
##   .. ..- attr(*, "names")= chr [1:4336366] "1" "2" "3" "4" ...

Calculos de Raster en R

Los calculos entre dos o mas rasters se realiza para crear una nueva salida. Por ejemplo si se esta interesado en mapear la altura de los arboles en toda una zona, se puede calcular la diferencia entre el modelo digital de superficie (DSM) y el modelo digital del terreno (DTM). como resultado el dataset se refiere al modelo de altura del canopi (CHM) y representa la altura actual de los arboles, edificios, etc, sin la influencia de la elevacion del suelo.

  • Ploteo Modelo Digital de Elevacion
ggplot() +
geom_raster(data = DEM_zona_df,
    aes(x = x, y = y,
    fill = DEM_zona)) +
scale_fill_gradientn(name = "Elevation", colors = terrain.colors(10))

  • Ploteo modelo digital de superficie

ggplot() + geom_raster(data = DEM_hill_df, aes(x = x, y = y, alpha = Hillshade_DEM1)) + scale_alpha(range = c(0.15, 0.65), guide = “none”)

ggplot() +
 geom_raster(data = Hill_df, aes(x = x, y = y, alpha = Hillshade_DEM1)) + 
    scale_alpha(range =  c(0.15, 0.65), guide = "none")

Dos formas de realizar calculos raster

Se puede calcular la diferencia entre dos rasters en dos formas diferentes:

  • Substrayendo los dos rasters en R utilizando raster math,

  • Ppara un proceso mas eficiente, particularmente si los rasters son grandes y el calculo se vuelve mas complejo, se utiliza la funcion overlay()

Matematica Raster y Modelos de altura del canopy

Es posible hacer calculos raster por sustracción, adicion, multiplicacion, etc entre dos rasters. En el mundo geoespacial, se le da el nombre de “raster math”

Resta del Modelo Digital del Terreno de el modelo digital de elevacion para crear el modelo de altura del canopy. Despues de la sustraccion, se crea un dataframe para plotear con ggplot.

CHM <- DEM - Hillshade

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

Plot del modelo de altura del canopy (CHM)

ggplot() +
      geom_raster(data = CHM_df , 
              aes(x = x, y = y, 
                   fill = layer)
              ) + 
     scale_fill_gradientn(name = "Altura Canopy", colors = terrain.colors(10))

Distribucion de los valores de la nueva capa creada de altura del canopy

ggplot(CHM_df) +
    geom_histogram(aes(layer))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

El rango esta entre 0 y 2600

Eficiencia de los calculos de raster: Funcion Overlay

Raster math, es una aproximacion apropiada para el calculo raster si:

  1. Los rasters que se esten utilizando utilizan un area pequeña en tamaño
  2. Los calculos que se realizan son simples

Sin embargo, raster math es poco eficiente en calculos mas complejos o archivos muy grandes. La funcion overlay() es mas eficiente cuando:

  1. Los rasters que se usan son grandes en tamaño.
  2. Los rasters estan almacenados como archivos individuales.
  3. Los calculos realizados son complejos.

La funcion overlay() toma dos o mas rasters y aplica una funcion utilizando metodos de procesamiento eficientes. La sintaxis es:

outputRaster <- overlay(raster1, raster2, fun=functionName)

Si los rasters estan almacenados como RasterStack o RasterBrick en R se debe usar calc(). Ya que overlay() no funciona con raster stack

Hora se va a realizar el mismo calculo de resta en que se utilizó raster math, pero ahora utilizando la funcion overlay().

CHM_ov <- overlay(DEM, Hillshade,
               fun = function(r1, r2) { return( r1 - r2) })
plot(CHM_ov)

Ahora se necesita convertir el nuevo objeto en un data frame para plotear con ggplot

CHM_ov_df <- raster::as.data.frame(CHM_ov, xy = TRUE)
CHM_ov_df <- na.omit(CHM_ov_df)
  • Plot del modelo de altura del canopy (CHM)
ggplot() +
      geom_raster(data = CHM_ov_df, 
              aes(x = x, y = y, 
                   fill = layer)
              ) + 
     scale_fill_gradientn(name = "Altura canopy", colors = terrain.colors(10))

Exportar un GeoTIFF

Para exportar el raster creado como un GeoTIFF se utiliza la funcion writeRaster().

Cuando se escribe un objeto raster a un archivo GeoTIFF, el nombre queda nombre-tiff. La funcion writeRaster() escribe por defecto el archivo de salida en su directorio de trabajo a menos que se especifique una ruta de archivo diferente.

Se especifica el formato de salida (“GTiff”), los Na (NAflag =-9999). Tambien se debe especificar a R sobre reescribir sobre un archivo con el mismo nombre (overwrite =TRUE)

writeRaster(CHM_ov, "CHM.tiff",
            format="GTiff",
            overwrite=TRUE,
            NAflag=-9999)