INTRODUCCIÓN

Qúe es R?

Nuevos avances en el campo de la informática han permitido el desarrollo de poderosas herramienta para la visualización, manejo y análisis de datos espaciales en áreas como la ecología del paisaje, ecología espacial y ciencias ambientales; cuya principal ventaja es la posibilidad de manejar en un solo “ambiente informático” múltiples tipos de datos, de distintos orígenes y prepararlos para su análisis, ya sean de carácter informativo, estadístico o más avanzados, por ejemplo, la predicción de nichos ecológicos de especies, sin necesidad de migrar los datos entre distintos tipos software.

El programa R, según su propia definición es “un lenguaje y un ambiente para la computación estadística y de gráficos”, está disponible en https://cran.r-project.org/ de forma gratuita y hace parte de las iniciativas de software de Acceso Abierto (en inglés, Open Access (OA)). R, como plataforma estadística en lenguaje de programación, se ha convertido en la última década en una de las herramientas más populares tanto para el uso de estadística clásica como para las técnicas de análisis más innovadoras en múltiples disciplinas.

Aunque R empezó como un lenguaje estadístico, ha evolucionado a un entorno de programación multipropósito que incluye funcionalidades de Sistemas de Información Geográfica y Ecología Espacial. R ofrece un número increíble de extensiones denominadas paquetes o librerías (Packages en inglés, ver https://cran.r-project.org/web/packages/index.html) que son desarrollados por la comunidad de usuarios de R.

Estos paquetes cubren todos los aspectos relacionados con la manipulación de datos, estadística, minería de datos o diseño de gráficos. Muchos de estos paquetes son creados con propósitos específicos, por ejemplo, para el análisis de índices de Biodiversidad o para la construcción de los modelos de nicho ecológico. Todos los paquetes disponibles en R han pasado por una rigurosa revisión para garantizar su aplicabilidad y correcto funcionamiento, estos esta agrupados por “grupos” de paquetes que tienen una utilidad en algún campo del conocimiento, los cuales están acompañados de comentarios y explicaciones, para nuestros objetivos puede revisar la siguiente vistas de tarea (https://cran.r-project.org/web/views/Spatial.html) que se han creado para reunir aquellos paquetes creados con fines de análisis espaciales.

Es importante recalcar que R tiene una comunidad muy activa, por lo que, haciendo las preguntas correctas rápidamente se puede encontrar la solución a los problemas que se presenten en el ámbito de la programación con R, lo que ha promovido que el número de usuarios de R en el área de las ciencias se incremente enormemente (Santana & Farfán, 2014).

INSTALANDO R Y RSTUDIO

R

En los siguientes enlaces puede encontrar los instaladores de R (tanto para Windows como para Mac):

R para Windows: https://cran.r-project.org/bin/windows/base/

R para Mac: https://cran.r-project.org/bin/macosx/

Siga el proceso de instalación de R una vez descargados los instaladores

RStudio

RStudio es una interface gráfica de R que, a su vez es el lenguaje de programación (Mas, 2013).

R Studio: https://rstudio.com/products/rstudio/download/#download

DATOS PARA LA PRÁCTICA

En el siguiente link podrá encontrar los datos necesarios para realizar la práctica:

https://drive.google.com/drive/folders/1inkcesNmIJt4p5eOuA2m9LIeJW_9GG8O?usp=sharing

Puede acceder a ellos en cualquier momento a través de su cuenta institucional (dominio @unal)

Dentro del archivo .rar descargado en el enlace de arriba se encuentra la siguiente información:

RSTUDIO

Aquí se presenta brevemente los diferentes paneles de RStudio que podemos observar al iniciar el programa por primera vez:

Interfaz de RStudio

Para ejecutar líneas de código dentro del editor de código seleccionamos las líneas deseadas y damos click en Run:

Ejecutar líneas de codigo

Dentro de esta guía podremos encontrar dos tipos de secciones que nos ayudará a identificar cuál es el código que ejecutaremos dentro del editor de código (panel con fondo en azul claro: 1) y cuál es la salida esperada en consola (panel con fondo gris: 2)

INSTALANDO Y CARGANDO PAQUETES NECESARIOS

Para instalar los paquetes necesarios para esta práctica ejecutamos el siguiente comando en la consola de RStudio (la instalación de paquetes es necesario hacerla una única vez.):

install.packages(“sf”,“spatialEco”,“dplyr”,“ggplot2”,“raster”,“rasterVis”)

Cuando se corren lineas de código que se demoran en ser ejecutadas aparecerá un hexagono rojo con la palabra STOP en color blanco (ver figura de abajo). Esto significa que se están ejecutando las intrucciones, sea paciente y espere, la instalación de paquetes podría ser demorada.

Para cargar los paquetes anteriormente descargados empleamos la función library(), junto con el nombre del paquete dentro de nuestro editor de código. Para crear un nuevo script en blanco sobre el que trabajar vamos a File > New File > R Script

library("sf")
library("spatialEco")
library("dplyr")
library("ggplot2")
library("raster")
library("rasterVis")

TRABAJANDO CON INFORMACIÓN EN FORMATO VECTORIAL

Una vez descargado el archivo .rar que contiene los datos necesarios para la práctica los extraemos en alguna ubicación de nuestra preferencia para así fijar un directorio de trabajo. Por ejemplo, en este caso los extraeremos dentro de una nueva carpeta llamada Datos practica que creamos en la unidad D. Por tanto, nuestro directorio de trabajo estará ubicado en la ruta: “D:/Datos practica”

A continuación, fijaremos nuestro directorio de trabajo dentro del entorno de R:

La función setwd() permite fijar nuestro directorio de trabajo, mientras qué getwd() sirve para consultar cual es el actual directorio

setwd("D:/Datos practica")

getwd()
## [1] "D:/Datos practica"

Para cargar el shapefile de departamentos de Colombia que está disponible en el directorio de trabajo usamos la función st_read() del paquete sf y guardamos el shapefile como objeto de R con el nombre de departamentos.colombia. Al cargar el shapefile, en consola nos aparecerá de forma automática información referente al archivo, en donde podremos ver que este consta por de 32 registros (features) y 2 campos o atributos (fields), entre otra información

departamentos.colombia <- st_read("Departamentos.shp")
## Reading layer `Departamentos' from data source `D:\Datos practica\Departamentos.shp' using driver `ESRI Shapefile'
## Simple feature collection with 32 features and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -9110557 ymin: -467986.2 xmax: -7443971 ymax: 1782881
## projected CRS:  WGS 84 / World Mercator

Con la función head() podemos examinar los primeros 3 registros de departamentos.colombia, que son Amazonas, Antioquia y Arauca, también podemos observar que el shapefile contiene dos atributos, el nombre del departamento y el área de este en km2

head(departamentos.colombia, 3)

De la misma forma cargamos el shapefile de incendios y lo guardamos como objeto de R con el nombre de incendios.2020, y examinamos los 3 primeros registros. Estos datos de incendios corresponden a anomalías térmicas detectadas por el sensor MODIS, disponibles para descarga gratuita en la base de datos FIRMS (Fire Information for Resource Management System) (https://firms.modaps.eosdis.nasa.gov/)

incendios.2020 <- st_read("incendios.shp")
## Reading layer `Incendios' from data source `D:\Datos practica\Incendios.shp' using driver `ESRI Shapefile'
## Simple feature collection with 9863 features and 3 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -8676865 ymin: -165836.4 xmax: -7502767 ymax: 1287783
## projected CRS:  WGS 84 / World Mercator
head(incendios.2020, 3)

Lo siguiente será graficar los shapefiles cargados. La forma más básica de graficar un objeto de R es usando la función plot(). Primeramente haremos el intento usando esta función con el archivo departamentos.colombia para ver qué pasa

plot(departamentos.colombia)

Automáticamente se generan dos gráficos, uno por cada atributo que contiene el shapefile de departamentos. Si deseamos graficar un unico atributo es necesario especificar cuál. A continuación, graficaremos únicamente el atributo de Area_km2. El argumento key.pos nos permite cambiar la localización de la barra de leyenda, y con el argumento axes = TRUE especificamos que aparezcan los valores de los ejes X y Y dentro del gráfico

plot(departamentos.colombia["Area_km2"], key.pos = 1, axes = TRUE)

Para graficar únicamente la geometría del objeto nos valimos de la función st_geometry(). A continuación, añadiremos en un único grafico los shapes de departamentos e incendios. Los argumentos col, pch y cex nos permiten cambiar el color, forma (puntos, triángulos, etc) y tamaño de la geometría de la capa de incendios, con el argumento ** add = TRUE** estamos indicando que se grafiquen los incendios encima de la capa de departamentos. Los argumentos main y axes=True sirven para agregar un título y ejes al gráfico.

plot(st_geometry(departamentos.colombia), main = "Incendios en Colombia (2020)", axes = TRUE)

plot(st_geometry(incendios.2020), col = "red", pch = 16, cex = 0.1, add = TRUE)

A continuación, realizaremos un análisis de carácter espacial. Queremos saber el número de incendios que ocurren dentro de los límites de cada departamento. Para esto usamos la función point.in.poly del paquete spatialEco. Esta función permite realizar la intersección entre polígonos y puntos, asignándole a cada punto los atributos del polígono que lo contiene, por tanto, cada punto del shapefile de incendios tendrá un nuevo atributo, en donde se indicará el departamento al que pertenece. Guardamos esta nueva capa con el nombre de incendios.departamentos

incendios.departamentos <- point.in.poly(departamentos.colombia, incendios.2020)

head(incendios.departamentos, 5)

Una vez que tenemos el nuevo conjunto de datos, donde a cada incendio se le asignó el atributo del nombre del departamento donde este se presentó es necesario construir una tabla con la frecuencia de incendios por departamento. Para este fin se usan las funciones table() que genera una tabla con el respectivo conteo por departamento y as.data.frame que convierte dicha tabla a un formato de filas y columnas para facilitar su manejo y visualización

departamentos.cuenta <- as.data.frame(table(incendios.departamentos$Nombre))

Ahora examinaremos los primeros 10 registros de la nueva tabla de conteo. Podemos ver que para cada departamento está asociado un valor de frecuencia, que es el número de incendios que se presentan en cada uno de ellos

head(departamentos.cuenta,10)

Sin embargo, solo deseamos contar con los 10 de departamentos que presentan la mayor frecuencia de fuegos. Para esto empleamos la función arrange() del paquete dplyr para reordenar la tabla en orden descendente de frecuencia. Lo siguiente es aplicar la función head() a nuestra tabla con el fin de seleccionar solo los primeros 10 registros

departamentos.cuenta <- arrange(departamentos.cuenta, -Freq)

departamentos.cuenta <- head(departamentos.cuenta, 10)

Lo siguiente será generar un gráfico de barras con nuestra última tabla de conteo de incendios por departamento. Para esto usaremos la función ggplot() del paquete ggplot2. Los argumentos de esta función serán data (la tabla a partir de la cual queremos realizar el gráfico) y aes para indicar cual columna queremos que este en el eje X y cual en el eje Y. A lo anterior debemos adicionar una línea con la función geom_bar() que determina el tipo de grafico que queremos, en este caso, de barras. Este grafico lo almacenaremos en un objeto llamado grafico.frecuencia1. Para visualizar el gráfico simplemente agregamos una línea con el nombre del objeto que lo contiene.

grafico.frecuencia1 <- ggplot(data=departamentos.cuenta, aes(x=Var1, y=Freq)) +  geom_bar(stat="identity")

grafico.frecuencia1

Cómo podemos se genera un gráfico muy básico. Para hacerlo más llamativo fácil y de interpretar vamos a generar un gráfico con más detalles. Para esto vamos a adicionar una serie de líneas y funciones. Dentro de la función geom_bar() podemos usar los argumentos color y fill para indicar de qué color deseamos que sean las barras, así como el color del borde de estas. Pruebe a usar diferentes combinaciones de colores.

grafico.frecuencia2 <- ggplot(data=departamentos.cuenta, aes(x=Var1, y=Freq)) + geom_bar(stat="identity", color="red", fill="black")

grafico.frecuencia2

también podemos usar la función labs() para agregar o cambiar a nuestro gráfico los títulos de los ejes X y Y, así como un título general para el gráfico

grafico.frecuencia3 <- ggplot(data=departamentos.cuenta, aes(x=Var1, y=Freq)) + geom_bar(stat="identity", color="red", fill="black") + labs(title="Los 10 departamentos con mayor numero de incendios en Colombia (2020)", x ="Departamentos", y = "Número de incendios")

grafico.frecuencia3

TRABAJANDO CON INFORMACIÓN EN FORMATO RASTER

En esta sección obtendremos las cifras de coberturas de la tierra presentes en un área protegida en particular. Primeramente, cargaremos el shape de áreas protegidas de Colombia. Este archivo lo podemos obtener desde el portal de la Base Mundial de Datos sobre Áreas Protegidas (WDPA : World Database on Protected Areas), disponibles en el portal https://www.protectedplanet.net

areas.protegidas <- st_read("AreasProtegidas.shp")
## Reading layer `AreasProtegidas' from data source `D:\Datos practica\AreasProtegidas.shp' using driver `ESRI Shapefile'
## Simple feature collection with 442 features and 2 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -9128128 ymin: -428910.6 xmax: -7534478 ymax: 1676941
## projected CRS:  WGS 84 / World Mercator
head(areas.protegidas, 10)

Podemos observar que el shapefile de áreas protegidas posee 2 campos o atributos de interés, el nombre del área (Nombre) y la categoría IUCN a la que esta pertenece (Cat_IUCN)

A continuación, graficaremos el shape de áreas protegidas para tener una idea de la distribución espacial de estas dentro del territorio colombiano

plot(st_geometry(areas.protegidas),  main = "Áreas protegidas de Colombia", axes = TRUE)

Lo siguiente será cargar el raster de coberturas y guardarlo como un objeto de R bajo el nombre de coberturas y graficarlo

coberturas <- raster('Coberturas.tif')

plot(coberturas)

La capa de coberturas de la tierra presenta valores entre 1 y 17, en dónde cada uno de ellos está asociado a un tipo de cobertura, de la siguiente forma:

Por el momento no nos preocuparemos por la visualización de esta capa, de eso nos ocuparemos más adelante. Lo siguiente será seleccionar un área protegida de interés dentro del shape de áreas protegidas de Colombia. En este caso nos interesaremos por el Parque Nacional Natural (PNN) Serranía de Chiribiquete.

Ya que conocemos el nombre del área protegida de nuestro de interés filtraremos el conjunto de datos a través del atributo Nombre, usando la función filter() del paquete dplyr

chiribiquete <- areas.protegidas %>% filter(Nombre == "Serranía de Chiribiquete")

plot(st_geometry(chiribiquete), main = "PNN Serranía de Chiribiquete", axes = TRUE)

Para cortar el raster de coberturas por los límites del PNN Chiribiquete son necesarios dos pasos:

mascara.chiribiquete <- mask(x = coberturas, mask = chiribiquete)

plot(mascara.chiribiquete, main = "Mascara")

corte.chiribiquete <- crop(x = mascara.chiribiquete, y = extent(chiribiquete))

plot(corte.chiribiquete, main = "Corte")

Ahora que tenemos las coberturas del PNN Chiribiquete vamos a calcular el área total de cada una de ellas dentro de los límites del parque. Para esto usamos la función freq() sobre el raster corte.chiribiquete y la función as.data.frame, lo que nos arrojara una tabla de frecuencia de pixeles para cada tipo de cobertura (los pixeles etiquetados con NA corresponden a pixeles vacíos, sin información)

coberturas.cuenta <- as.data.frame(freq(corte.chiribiquete))

coberturas.cuenta

A esta nueva tabla podemos agregarle un campo para calcular el área en Km2, para esto hay que tomar en consideración que la resolución espacial del raster de coberturas es de 500m (es decir, cada pixel representa un área de 250.000m, o 0,25km2), por tanto debemos multiplicar el número de pixeles por 0,25 para obtener el área total. Para facilitar la interpretación de la tabla podemos asignar el texto que corresponde a cada valor de tipo de coberturas

coberturas.cuenta$Area_km2 <- coberturas.cuenta$count * 0.25

coberturas.cuenta$Cobertura <- c("Bosques perennes de hoja ancha","Bosques caducifolios de hoja ancha","Matorrales cerrados","Sabanas arboladas","Sabanas","Pastizales","Humedales permanentes","Sin información")

coberturas.cuenta

En el paso anterior nos dimos cuenta de que dentro de los límites del PNN tenemos 7 diferentes tipos de coberturas (valores 2,4,6,8,9,10,11), que nos servirá para mejorar el gráfico de coberturas para el PNN. El siguiente código permite graficar de forma más adecuada las coberturas:

corte.chiribiquete.a <- ratify(corte.chiribiquete)

corte.chiribiquete.b <- levels(corte.chiribiquete.a)[[1]]

corte.chiribiquete.b$legend <- c("Bosques perennes de hoja ancha","Bosques caducifolios de hoja ancha","Matorrales cerrados","Sabanas arboladas","Sabanas","Pastizales","Humedales permanentes")

levels(corte.chiribiquete.a) <- corte.chiribiquete.b

levelplot(corte.chiribiquete.a, col.regions = c("green","red","black","yellow","white","brown","blue"))

REFERENCIAS