Recursos de interés

Preparativos

Paquetes

# Paquete para manejo de datos vectoriales
library(sf)

# Paquete de Tidyverse para manipulación de datos
library(dplyr)

# Paquete para mapas en la Web
library(leaflet)

Datos para ejemplos

# Registros de presencia de Ara ambiguus (lapa verde) en Costa Rica
cr_ara_ambiguus <-
  st_read(
    "https://raw.githubusercontent.com/geoprocesamiento-2020i/datos/master/biodiversidad/registros-presencia/cr/cr-ara-ambiguus.geojson"
  )

El modelo de datos raster

El modelo de datos raster usualmente consiste de un encabezado y de una matriz con celdas (también llamadas pixeles) de un mismo tamaño. El encabezado define el sistema de referencia de coordenadas (CRS), la extensión y el punto de origen de una capa raster. Por lo general, el origen se ubica en la esquina inferior izquierda o en la esquina superior izquierda de la matriz. La extensión se define mediante el número de filas, el número de columnas y el tamaño (resolución) de la celda.

Cada celda tiene una identificación (ID) y almacena un único valor, el cual puede ser numérico o categórico, como se muestra en la figura 1.

Figura 1: El modelo raster: (A) ID de las celdas, (B) valores de las celdas, (C) mapa raster de colores. Imagen de Robin Lovelace et al. (https://geocompr.robinlovelace.net/spatial-class.html#raster-data)

A diferencia del modelo vectorial, el modelo raster no necesita almacenar todas las coordenadas de cada geometría (i.e. las esquinas de las celdas), debido a que la ubicación de cada celda puede calcularse a partir de la información contenida en el encabezado. Esta simplicidad, en conjunto con el álgebra de mapas, permiten que el procesamiento de datos raster sea mucho más eficiente que el procesamiento de datos vectoriales. Por otra parte, el modelo vectorial es mucho más flexible en cuanto a las posibilidades de representación de geometrías y almacenamiento de valores, por medio de múltiples elementos de datos.

Los mapas raster generalmente almacenan fenómenos continuos como elevación, precipitación, temperatura, densidad de población y datos espectrales. También es posible representar mediante raster datos discretos, tales como tipos de suelo o clases de cobertura de la tierra, como se muestra en la figura 2.

Figura 2: Ejemplos de mapas raster continuos y categóricos. Imagen de Robin Lovelace et al. (https://geocompr.robinlovelace.net/spatial-class.html#raster-data)

El paquete raster

El paquete raster proporciona funciones para la lectura, escritura, manipulación, análisis y modelado de datos raster.

El paquete rgdal provee enlaces a las bibliotecas GDAL y PROJ.

Instalación, carga y documentación

Instalación

# Instalación del paquete raster
install.packages("raster")

# Instalación del paquete raster
install.packages("rgdal")

Carga

# Carga del paquete raster
library(raster)

# Carga del paquete rgdal
library(rgdal)

Ayuda

# Ayuda sobre paquete raster
help("raster-package")

Documentación:
- The raster package – R Spatial
- raster: Geographic Data Analysis and Modeling

Formatos raster

El paquete raster puede leer y escribir en una gran cantidad de formatos raster a través de la biblioteca GDAL.

Uso básico

Para ejemplificar el uso del paquete raster, se accederá a la base de datos climáticos WorldClim, mediante la función getData().

# Directorio de trabajo (DEBE USARSE UN DIRECTORIO EXISTENTE EN EL DISCO)
setwd("c:/users/mfvargas/")

# Consulta del directorio de trabajo
getwd()

## [1] "c:/users/mfvargas"

# Datos de altitud
altitude <- getData("worldclim", var="alt", res=.5, lon=-84, lat=10)

# Datos de altitud recortados para los límites aproximados de Costa Rica
cr_altitude <- crop(altitude, extent(-86, -82.3, 8, 11.3))  

# Resumen de información básica de la capa raster
cr_altitude

## class      : RasterLayer 
## dimensions : 396, 444, 175824  (nrow, ncol, ncell)
## resolution : 0.008333333, 0.008333333  (x, y)
## extent     : -86, -82.3, 8, 11.3  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : alt_23 
## values     : -21, 3732  (min, max)

Tipos de datos

El paquete raster define dos tipos de datos principales:

RasterLayer: almacena una sola capa.

# Clase de cr_altitude
class(cr_altitude)

## [1] "RasterLayer"
## attr(,"package")
## [1] "raster"

# Cantidad de capas
nlayers(cr_altitude)

## [1] 1

# Mapeo
plot(cr_altitude)

RasterBrick: almacena varias capas, las cuales generalmente provienen de un mismo origen (ej. una imagen satelital).

# Datos de precipitación
precipitation <- getData("worldclim", var="prec", res=.5, lon=-84, lat=10)

# Datos de precipitación recortados para los límites aproximados de Costa Rica
cr_precipitation <- crop(precipitation, extent(-86, -82.3, 8, 11.3))

# Clase de cr_precipitacion
class(cr_precipitation)

## [1] "RasterBrick"
## attr(,"package")
## [1] "raster"

# Cantidad de capas
nlayers(cr_precipitation)

## [1] 12

# Lista e información sobre las capas
unlist(cr_precipitation)

## class      : RasterBrick 
## dimensions : 396, 444, 175824, 12  (nrow, ncol, ncell, nlayers)
## resolution : 0.008333333, 0.008333333  (x, y)
## extent     : -86, -82.3, 8, 11.3  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : prec1_23, prec2_23, prec3_23, prec4_23, prec5_23, prec6_23, prec7_23, prec8_23, prec9_23, prec10_23, prec11_23, prec12_23 
## min values :        0,        0,        0,        9,      104,      179,      122,      154,       90,       140,        79,         7 
## max values :      500,      293,      283,      372,      618,      659,      807,      673,      790,       829,       801,       752

# Mapeo
plot(cr_precipitation)

RasterStack: también almacena varias capas, pero a diferencia de RasterBrick, estas provienen de diferentes fuentes.

Confección de mapas

Ejemplo de mapa en Leaflet

# Paleta de colores
pal <- colorNumeric(
  c("#0C2C84", "#41B6C4", "#FFFFCC"), 
  values(cr_altitude), 
  na.color = "transparent"
)

# Mapa web
m <- leaflet() %>% 
  addProviderTiles(providers$Esri.WorldImagery, group = "Imágenes de ESRI") %>% 
  addProviderTiles(providers$Stamen.TonerLite, group = "Stamen Toner Lite") %>%
  addProviderTiles(providers$OpenStreetMap.Mapnik, group = "OpenStreetMap") %>%
  addCircleMarkers(data = cr_ara_ambiguus,
                   stroke = F,
                   radius = 4,
                   fillColor = 'green',
                   fillOpacity = 1,
                   group = "Ara ambiguus",
                   popup = paste(cr_ara_ambiguus$locality, 
                                 cr_ara_ambiguus$year, 
                                 sep = '<br/>'
                           )
  ) %>%
  addRasterImage(cr_altitude, 
                 colors = pal, 
                 opacity = 0.8, 
                 group = "Altitud"
  ) %>%
  addLayersControl(
    baseGroups = c("OpenStreetMap", "Stamen Toner Lite", "Imágenes de ESRI"),
    overlayGroups = c("Altitud", "Ara ambiguus"),
    options = layersControlOptions(collapsed = FALSE)    
  ) %>%
  addLegend(pal = pal, 
            values = values(cr_altitude), 
            title = "Altitud"
  ) %>%
  addMiniMap(
    toggleDisplay = TRUE,
    position = "bottomleft",
    tiles = providers$Stamen.TonerLite
  )

# Despliegue del mapa
m

Análisis

# Directorio de trabajo (DEBE USARSE UN DIRECTORIO EXISTENTE EN EL DISCO)
setwd("c:/users/mfvargas/")

# Consulta del directorio de trabajo
getwd()

# Datos de altitud
altitud <- getData("worldclim", var="alt", res=.5, lon=-84, lat=10)

# Datos de altitud recortados para los límites aproximados de Costa Rica
cr_altitud_aproximada <- 
  altitud %>%
  crop(extent(-86, -82.3, 8, 11.3))
plot(cr_altitud_aproximada)

# Recorte para el contorno de Costa Rica
# 1. Capa de provincias
url_base_wfs_ign_5mil <- "http://geos.snitcr.go.cr/be/IGN_5/wfs?"
solicitud_provincias_wfs <- "request=GetFeature&service=WFS&version=2.0.0&typeName=IGN_5:limiteprovincial_5k&outputFormat=application/json"
cr_provincias <-
  st_read(paste0(url_base_wfs_ign_5mil, solicitud_provincias_wfs)) %>%
  st_simplify(dTolerance = 1000) %>%
  st_transform(4326)

# 2. Recorte
cr_altitud <-
  altitud %>%
  crop(cr_provincias) %>%
  mask(cr_provincias)

# 3. Mapeo de todo el territorio
plot(cr_altitud)

# 4. Mapeo de solo la parte continental
plot(cr_altitud, ext=extent(-86, -82.3, 8, 11.3))

# Registros de presencia de Ara ambiguus (lapa verde) en Costa Rica
cr_ara_ambiguus <-
  st_read(
    "https://raw.githubusercontent.com/geoprocesamiento-2020i/datos/master/biodiversidad/registros-presencia/cr/cr-ara-ambiguus.geojson"
  )

# Registros de presencia de Pharomachrus mocinno (quetzal) en Costa Rica
cr_pharomachrus_mocinno <-
  st_read(
    "https://raw.githubusercontent.com/geoprocesamiento-2020i/datos/master/biodiversidad/registros-presencia/cr/cr_pharomachrus_mocinno.geojson"
  )

# Mapeo
plot(cr_altitud, ext=extent(-86, -82.3, 8, 11.3), axes=TRUE, graticule=TRUE, reset=FALSE)
plot(cr_ara_ambiguus, col="green", add=TRUE)
plot(cr_pharomachrus_mocinno, col="black", add=TRUE)

# Mapa de altitud del área de Ara ambiguus
plot(cr_altitud, 
     ylim=c(min(cr_ara_ambiguus$decimalLatitude), 
            max(cr_ara_ambiguus$decimalLatitude)
          ), 
     xlim=c(min(cr_ara_ambiguus$decimalLongitude), 
            max(cr_ara_ambiguus$decimalLongitude)
          ),
     reset=FALSE
)
plot(cr_ara_ambiguus, col="black", add=TRUE)

# Mapa de altitud del área de Pharomachrus mocinno
plot(cr_altitud, 
     ylim=c(min(cr_pharomachrus_mocinno$decimalLatitude), 
            max(cr_pharomachrus_mocinno$decimalLatitude)
          ), 
     xlim=c(min(cr_pharomachrus_mocinno$decimalLongitude), 
            max(cr_pharomachrus_mocinno$decimalLongitude)
          ),
     reset=FALSE
)
plot(cr_pharomachrus_mocinno, col="black", add=TRUE)

# Altitudes de los puntos de Ara ambiguus
alt_ara_ambiguus <- extract(cr_altitud, cr_ara_ambiguus)
# Impresión 
alt_ara_ambiguus

# Altitudes de los puntos de Pharomachrus mocinno
alt_pharomachrus_mocinno <- extract(cr_altitud, cr_pharomachrus_mocinno)
# Impresión 
alt_pharomachrus_mocinno

Ejercicios

Genere un mapa de altitud para la provincia de Heredia y otro para la provincia de San José.
Genere un mapa de altitud para una especie diferente a las utilizadas durante la lección. Puede utilizar datos de la Infraestructura Mundial de Información en Biodiversidad (GBIF).

El lenguaje de programación R - Datos raster

Manuel Vargas

2020-06-11