En este ejercicio se muestra cómo se grafíca un archivo raster con diferentes reglas de color usando ggplot2.

Carga de librerías necesarias:

library(raster)
## 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.1.2, released 2016/10/24
##  Path to GDAL shared files: /usr/share/gdal/2.1
##  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: (autodetected)
##  Linking to sp version: 1.2-7
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
## 
##     intersect, select, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
GDALinfo("/home/jorge/Documentos/Maestria/PRA/Dataset_GEE/DEM_Bosavita_10m.tif")
## rows        1808 
## columns     1732 
## bands       1 
## lower left origin.x        -73.56242 
## lower left origin.y        5.256726 
## res.x       1.078094e-05 
## res.y       1.07873e-05 
## ysign       -1 
## oblique.x   0 
## oblique.y   0 
## driver      GTiff 
## projection  +proj=longlat +datum=WGS84 +no_defs 
## file        /home/jorge/Documentos/Maestria/PRA/Dataset_GEE/DEM_Bosavita_10m.tif 
## apparent band summary:
##    GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64           TRUE           0          1       1732
## apparent band statistics:
##       Bmin     Bmax   Bmean      Bsd
## 1 2760.348 3045.331 2914.69 50.75652
## Metadata:
## AREA_OR_POINT=Area

Abrir el archivo raster en R

DSM_BVTA <- 
  raster("/home/jorge/Documentos/Maestria/PRA/Dataset_GEE/DEM_Bosavita_10m.tif")

DSM_BVTA
## class       : RasterLayer 
## dimensions  : 1808, 1732, 3131456  (nrow, ncol, ncell)
## resolution  : 1.078094e-05, 1.07873e-05  (x, y)
## extent      : -73.56242, -73.54375, 5.256726, 5.27623  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /home/jorge/Documentos/Maestria/PRA/Dataset_GEE/DEM_Bosavita_10m.tif 
## names       : DEM_Bosavita_10m 
## values      : 2760.348, 3045.331  (min, max)
DSM_BVTA_df <- raster::as.data.frame(DSM_BVTA, xy = TRUE)

Ahora se puede ver la imagen como un dataframe:

str(DSM_BVTA_df)
## 'data.frame':    3131456 obs. of  3 variables:
##  $ x               : num  -73.6 -73.6 -73.6 -73.6 -73.6 ...
##  $ y               : num  5.28 5.28 5.28 5.28 5.28 ...
##  $ DEM_Bosavita_10m: num  NA NA NA NA NA NA NA NA NA NA ...

Graficar el archivo raster con rangos de color:

En el ejercicio 1 se graficó el archivo raster utilizando una rampa continua de colores. En este ejercicio usaremos reglas o rangos de colores para representar nuestro DSM.

library(ggplot2)

DSM_BVTA_df <- DSM_BVTA_df %>%
                mutate(fct_elevation = cut(DEM_Bosavita_10m, breaks = 3))

ggplot() +
    geom_bar(data = DSM_BVTA_df, aes(fct_elevation))

A continuación podemos ver los valores de corte de nuestros rangos:

unique(DSM_BVTA_df$fct_elevation)
## [1] <NA>                (2.76e+03,2.85e+03] (2.85e+03,2.95e+03]
## [4] (2.95e+03,3.05e+03]
## Levels: (2.76e+03,2.85e+03] (2.85e+03,2.95e+03] (2.95e+03,3.05e+03]

Utilizando las funciones group_by() y summarize() de la librería dplyr podemos ver la cantidad de datos en cada rango:

DSM_BVTA_df %>%
        group_by(fct_elevation) %>%
        summarize(counts = n())

Ahora podemos graficar nuestro raster:

ggplot() +
 geom_raster(data = DSM_BVTA_df , aes(x = x, y = y, fill = fct_elevation))

También podemos configurar los rangos manualmente:

custom_bins <- c(2700, 2800, 2900, 3100)

DSM_BVTA_df <- DSM_BVTA_df %>%
                mutate(fct_elevation_2 = cut(DEM_Bosavita_10m, breaks = custom_bins))

unique(DSM_BVTA_df$fct_elevation_2)
## [1] <NA>              (2.8e+03,2.9e+03] (2.7e+03,2.8e+03] (2.9e+03,3.1e+03]
## Levels: (2.7e+03,2.8e+03] (2.8e+03,2.9e+03] (2.9e+03,3.1e+03]

Ahora podemos graficar nuevamente nuestros rangos:

ggplot() +
    geom_bar(data = DSM_BVTA_df, aes(fct_elevation_2))

Y podemos revisar nuevamente la cantidad de valores en cada rango:

DSM_BVTA_df %>%
        group_by(fct_elevation_2) %>%
        summarize(counts = n())

Con estos rangos podemos graficar nuestro raster:

ggplot() +
 geom_raster(data = DSM_BVTA_df , aes(x = x, y = y, fill = fct_elevation_2))

Usando la función terrain.colors() podemos especificar colores más apropiados para una representación de topografía y mejorar la visualización de nuestro raster:

ggplot() +
 geom_raster(data = DSM_BVTA_df , aes(x = x, y = y,
                                      fill = fct_elevation_2)) + 
    scale_fill_manual(values = terrain.colors(3))
## Warning: Removed 1545995 rows containing missing values (geom_raster).

Opciones de formato del gráfico

También podemos agregar algunas opciones de formato al gráfico:

my_col <- terrain.colors(3)

ggplot() +
 geom_raster(data = DSM_BVTA_df , aes(x = x, y = y,
                                      fill = fct_elevation_2)) + 
    scale_fill_manual(values = my_col, name = "Elevación")
## Warning: Removed 1545995 rows containing missing values (geom_raster).

Se pueden eliminar las etiquetas de los ejes:

ggplot() +
 geom_raster(data = DSM_BVTA_df , aes(x = x, y = y,
                                      fill = fct_elevation_2)) + 
    scale_fill_manual(values = my_col, name = "Elevación") +
    theme(axis.title.x = element_blank(),
          axis.title.y = element_blank())
## Warning: Removed 1545995 rows containing missing values (geom_raster).

Visualización del DSM con efecto 3D

Primero se debe cargar un mapa de sombras del DSM:

DSM_hill_BVTA <-
  raster("/home/jorge/Documentos/Maestria/PRA/Dataset_GEE/HILL_Bosavita_10m.tif")

DSM_hill_BVTA
## class       : RasterLayer 
## dimensions  : 1808, 1732, 3131456  (nrow, ncol, ncell)
## resolution  : 1.078094e-05, 1.07873e-05  (x, y)
## extent      : -73.56242, -73.54375, 5.256726, 5.27623  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /home/jorge/Documentos/Maestria/PRA/Dataset_GEE/HILL_Bosavita_10m.tif 
## names       : HILL_Bosavita_10m 
## values      : 1, 182  (min, max)

Ahora convertimos el mapa de sombras a dataframe:

DSM_hill_BVTA_df <- raster::as.data.frame(DSM_hill_BVTA, xy = TRUE) 

str(DSM_hill_BVTA_df)
## 'data.frame':    3131456 obs. of  3 variables:
##  $ x                : num  -73.6 -73.6 -73.6 -73.6 -73.6 ...
##  $ y                : num  5.28 5.28 5.28 5.28 5.28 ...
##  $ HILL_Bosavita_10m: int  NA NA NA NA NA NA NA NA NA NA ...

Ahora graficamos el mapa de sombras:

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

Finalmente podemos graficar nueestro DSM con un efecto 3D:

ggplot() +
    geom_raster(data = DSM_BVTA_df , 
                aes(x = x, y = y, 
                     fill = DEM_Bosavita_10m)) + 
  
    geom_raster(data = DSM_hill_BVTA_df, 
                aes(x = x, y = y, 
                  alpha = HILL_Bosavita_10m)) +
  
    scale_fill_viridis_c() +
  
    scale_alpha(range = c(0.15, 0.65), guide = "none") +
  
    ggtitle("DSM con mapa de sombras") +
    coord_equal()