- Ploteo de Datos Raster en R
Ploteo de datos Raster 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 parte del tutorial de Introduccion de Datos Geoespaciales con R
Introducción
La imagen que se va a trabajar es la misma imagen utilizada en el primer tutorial que pertenece a una zona en el municipio de Villavicencio Meta.
Cargar librerias
library(raster)## Warning: package 'raster' was built under R version 3.4.4## Loading required package: splibrary(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-7library(dplyr)## Warning: package 'dplyr' was built under R version 3.4.4##
## 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, unionlibrary(ggplot2)Ubicacion del directorio y cargar la imagen de trabajo
# Seleccion de directorio de trabajo
setwd("C:\\Maestria en Geomatica\\Percepcion remota avanzada\\Talleres")
# Cargar imagen
info_B2 <- raster("Landsat_recorte_B2.TIF")Ploteo de Datos utilizando espacios
El ploteo de los datos se puede realizar con simbolos o colores de acuerdo con el rango de los valores para mejorar la claridad y visibilidadd del ploteo. Para esto se necesita decirle a ggplot en cuantos grupos de van a dividir los datos, y donde deben estar los descansos. Para tomar estas decisiones es util explorar primero la distribucion de los datos mediante diagrama de barras. Para esto se utiliza la funcion dplyr´s mulate() combinada con cut() para dividir los datos en 3 contenedores.
La funcion na.omith() sirve para eliminar los valores de NA de la imagen
info_B2_df <- raster::as.data.frame(info_B2, xy = TRUE)
info_B2_df <- na.omit(info_B2_df)
str(info_B2_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" ...info_B2_df <- info_B2_df %>%
mutate(ND = cut(Landsat_recorte_B2, breaks = 3))## Warning: package 'bindrcpp' was built under R version 3.4.4ggplot() +
geom_bar(data = info_B2_df, aes(ND))
Si se quiere conocer los valores de corte para los grupos, se puede solicitar los valores unicos de altura
unique(info_B2_df$ND)## [1] (-37.6,1.25e+04] (1.25e+04,2.51e+04] (2.51e+04,3.77e+04]
## Levels: (-37.6,1.25e+04] (1.25e+04,2.51e+04] (2.51e+04,3.77e+04]Para contar los valores en cada grupo se utilizan las funciones dplyr, group_by() y summarize().
info_B2_df %>%
group_by(info_B2_df$ND) %>%
summarize(counts = n())Para personalizar los valores de corte de cada grupo se debe utilizar la funcion mutate() a la cual se le da un vector numerico de puntos de ruptura en vez de darle la cantidad de rupturas que se quieren, esto para redondear los valores de corte de los grupos.
custom_bins <- c( 1000, 5000, 10000, 37000)
info_B2_df <- na.omit(info_B2_df)
summary(info_B2_df)## x y Landsat_recorte_B2
## Min. :616800 Min. :397530 Min. : 0
## 1st Qu.:645000 1st Qu.:421770 1st Qu.: 8647
## Median :668730 Median :433470 Median : 8922
## Mean :668774 Mean :434112 Mean : 8548
## 3rd Qu.:691800 3rd Qu.:446520 3rd Qu.: 9197
## Max. :723870 Max. :476040 Max. :37616
## ND
## (-37.6,1.25e+04] :4986743
## (1.25e+04,2.51e+04]: 24380
## (2.51e+04,3.77e+04]: 1389
##
##
## str(info_B2_df)## 'data.frame': 5012512 obs. of 4 variables:
## $ x : num 635820 635850 635880 635910 635940 ...
## $ y : num 476040 476040 476040 476040 476040 ...
## $ Landsat_recorte_B2: num 7796 7920 7951 8034 8240 ...
## $ ND : Factor w/ 3 levels "(-37.6,1.25e+04]",..: 1 1 1 1 1 1 1 1 1 1 ...info_B2_df <- (info_B2_df %>%
mutate(ND_2 = cut(Landsat_recorte_B2, breaks = custom_bins)))
unique(info_B2_df$ND_2)## [1] (5e+03,1e+04] (1e+04,3.7e+04] <NA>
## Levels: (1e+03,5e+03] (5e+03,1e+04] (1e+04,3.7e+04]Ploteo del grafico de barras utilizando nuevos grupos
ggplot() +
geom_bar(data = info_B2_df, aes(ND_2))
Tambien se puede obtener el numero de valores en cada grupo de manera similar que el anterior
info_B2_df %>%
group_by(ND_2) %>%
summarize(counts = n())Se pueden utilizar esos grupos para plotear los datos raster personalizado, utilizando un color difernte para cada grupo.
ggplot() +
geom_raster(data = info_B2_df , aes(x = x, y = y, fill = ND_2))
El ploteo utiliza los colores predeterminados por ggplot para objetos raster. Se puede especificar un colores personalizados para obtener un ploteo mas agradable. En R existe un conjunto de colores para ploteo en la funcion terrain.colors(). Como se tienen 3 contenedores, se pueden utilizar los primeros tres colores del terreno.
terrain.colors(3)## [1] "#00A600FF" "#ECB176FF" "#F2F2F2FF"La funcion Terrarin.colors() son una cadena de caracteres que representan un color. Para utilizarlos se pude utilizar la funcion scale_fill_manual()
ggplot() +
geom_raster(data = info_B2_df , aes(x = x, y = y,
fill = ND_2)) +
scale_fill_manual(values = terrain.colors(3))## Warning: Removed 285036 rows containing missing values (geom_raster).
Otros formatos para plotear
Si se necesitan crear multiples ploteos utilizando el mismo color de la paleta, se puece crear un objeto en R (my_col) para un conjunto de colores que deseo utilizar. SE puede cambiar rapidamente la paleta de colores a traves de todos los plots para modificar el objeto my_col, en cada ploteo individual.
Se puede etiquetar los ejes x y y del plot utilizando xlab y ylab. Tambien se puede dar a la leyenda con un titulo mas significativo dando un valor al argumento del nombre en la funcion scale_fill_manual().
my_col <- terrain.colors(3)
ggplot() +
geom_raster(data = info_B2_df , aes(x = x, y = y,
fill = ND_2)) +
scale_fill_manual(values = my_col, name = "Numeros Digitales")## Warning: Removed 285036 rows containing missing values (geom_raster).
Tambien se pueden quitar los ejes utilizando element_blank()en la parte relevante de la funcion theme().
ggplot() +
geom_raster(data = info_B2_df , aes(x = x, y = y,
fill = ND_2)) +
scale_fill_manual(values = my_col, name = "Numeros Digitales") +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank())## Warning: Removed 285036 rows containing missing values (geom_raster).
Capas de un raster
Se pueden aplicar capas a un raster en la parte superior con un raster con la misma area y utilizar transparencia para crear un efecto de sombra en 3 dimensiones. Un hillshade es un raster que mapea las sombras y texturas que se pueden ver sobre un terreno. Se puede presonalizar el color haciendo un ploteo en escala de grises
DEM_hill <- raster("C:\\Maestria en Geomatica\\Percepcion remota avanzada\\Talleres\\Hillshade_DEM1.tif")
DEM_hill## class : RasterLayer
## dimensions : 988, 1081, 1068028 (nrow, ncol, ncell)
## resolution : 91.12718, 91.12718 (x, y)
## extent : 1022983, 1121491, 889190.4, 979224 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## data source : C:\Maestria en Geomatica\Percepcion remota avanzada\Talleres\Hillshade_DEM1.tif
## names : Hillshade_DEM1
## values : 0, 254 (min, max)Se convierte el nuevo raster a dataframe y se eliminan los NA
DEM_hill_df <- raster::as.data.frame(DEM_hill, xy = TRUE)
DEM_hill_df <- na.omit(DEM_hill_df)
str(DEM_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" ...summary(DEM_hill_df)
Ploteo del hillshade, alpha hace referencia a la escala de transparencia
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")
Se puede sobreponer otro raster sobre el hillshade utilizando otra funcion llamada geom_raster().
ggplot() +
geom_raster(data = info_B2_df ,
aes(x = x, y = y,
fill = Landsat_recorte_B2)) +
geom_raster(data = DEM_hill_df,
aes(x = x, y = y,
alpha = Hillshade_DEM1)) +
scale_fill_gradient() +
scale_alpha(range = c(0.15, 0.65), guide = "none") +
ggtitle("Numeros digitales / elevacion del terreno ") +
coord_equal()
scale_fill_manual(values = my_col, name = "Numeros Digitales") +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank())## NULL