13 Introducción al manejo de datos geoespaciales en R

13.1 Resumen

El lenguaje de programación R cuenta con varios paquetes para manejo de datos geoespaciales, entre los que están el paquete sf para datos vectoriales, el paquete terra para datos raster y el paquete leaflet para mapas en la Web.

13.2 Trabajo previo

13.2.1 Lecturas

Lovelace, R., Nowosad, J., & Münchow, J. (2019). Geocomputation with R (capítulos 1 y 2). CRC Press. https://geocompr.robinlovelace.net/

13.3 Preparativos

13.3.1 Carga de paquetes

# Paquete para transformación de datos
library(dplyr)

13.4 Introducción

La comunidad de programadores de R ha desarrollado un conjunto de paquetes para el manejo de datos geoespaciales, tanto en formatos vectoriales como raster. Algunos de los principales de estos paquetes son:

  • El paquete sf. Ofrece un conjunto de funciones para el manejo de datos vectoriales, de acuerdo con el estándar Simple Features.

  • El paquete terra. Implementa un conjunto de funciones para el manejo de datos raster. Es una reimplementación del paquete raster.

  • El paquete leaflet. Es una implementación en R de la biblioteca Leaflet para el lenguaje JavaScript para la programación de mapas interactivos en páginas web.

Algunos paquetes de graficación estadística, como ggplot2 y plotly, también cuentan con capacidades para visualización de datos geoespaciales.

En CRAN Task View: Analysis of Spatial Data, puede encontrarse un resumen detallado de los paquetes geoespaciales de R.

13.5 Datos vectoriales

13.5.1 El modelo vectorial

El modelo vectorial de datos está basado en puntos localizados en un sistema de referencia de coordenadas (CRS). Los puntos individuales pueden representar objetos independientes (ej. postes eléctricos, cabinas telefónicas) o pueden también agruparse para formar geometrías más complejas como líneas (ej. ríos, caminos) o polígonos (ej. fincas, países, provincias). Por lo general, los puntos tienen solo dos dimensiones (x, y), a las que se les puede agregar una tercera dimensión z, usualmente correspondiente a la altitud sobre el nivel del mar.

13.5.2 El estándar Simple Features

Simple Features (o Simple Feature Access) es un estándar abierto de la Organización Internacional de Estandarización (ISO) y del Open Geospatial Consortium (OGC) que especifica un modelo común de almacenamiento y acceso para geometrías de dos dimensiones (líneas, polígonos, multilíneas, multipolígonos, etc.). El estándar es implementado por muchas bibliotecas y bases de datos geoespaciales como sf, Fiona, GDAL, PostgreSQL/PostGIS, SQLite/SpatiaLite, Oracle Spatial y Microsoft SQL Server, entre muchas otras.

La especificación define 17 tipos de geometrías, de las cuales siete son las más comúnmente utilizadas. Estas últimas se muestran en la figura 13.1.

Tipos de geometrías de Simple Features más usadas. Imagen de [Robin Lovelace et al.](https://geocompr.robinlovelace.net/spatial-class.html#vector-data)

Figure 13.1: Tipos de geometrías de Simple Features más usadas. Imagen de Robin Lovelace et al.

13.5.3 El paquete sf

El paquete sf (de Simple Features) de R implementa los modelos de datos de las geometrías de tipo vectorial: puntos, líneas, polígonos, sus versiones múltiples y las colecciones de geometrías. Está basado en bibliotecas de sofware ampliamente utilizadas en aplicaciones geoespaciales:

sf provee acceso, desde un mismo paquete de R, a la funcionalidad de estas tres bibliotecas, proporcionando así una interfaz unificada para leer y escribir datos geoespaciales mediante GDAL, realizar operaciones con geometrías mediante GEOS y efectuar transformaciones entre sistemas de coordenadas mediante PROJ.

En sf, los conjuntos de datos geoespaciales se almacenan en objetos de una clase también llamada sf, los cuales son data frames que contiene una columna especial para las geometrías. Esta columna se denomina generalmente geom o geometry (aunque pueden tener cualquier otro nombre). El manejo de datos geoespaciales como data frames permite manipularlos con las funciones ya desarrolladas para este tipo de datos y con la misma forma de referenciar las filas (observaciones) y las columnas (variables).

13.5.3.1 Instalación y carga

# Instalación de sf
install.packages("sf")
# Carga de sf
library(sf)

13.5.3.2 Métodos

La lista de métodos (i.e. funciones) de la clase sf puede obtenerse a través de la función methods():

# Métodos de la clase sf
methods(class = "sf")
#>  [1] [                     [[<-                 
#>  [3] $<-                   aggregate            
#>  [5] anti_join             arrange              
#>  [7] as.data.frame         cbind                
#>  [9] coerce                dbDataType           
#> [11] dbWriteTable          distinct             
#> [13] dplyr_reconstruct     filter               
#> [15] full_join             group_by             
#> [17] group_split           identify             
#> [19] initialize            inner_join           
#> [21] left_join             merge                
#> [23] mutate                plot                 
#> [25] print                 rbind                
#> [27] rename                right_join           
#> [29] rowwise               sample_frac          
#> [31] sample_n              select               
#> [33] semi_join             show                 
#> [35] slice                 slotsFromS3          
#> [37] st_agr                st_agr<-             
#> [39] st_area               st_as_s2             
#> [41] st_as_sf              st_as_sfc            
#> [43] st_bbox               st_boundary          
#> [45] st_buffer             st_cast              
#> [47] st_centroid           st_collection_extract
#> [49] st_convex_hull        st_coordinates       
#> [51] st_crop               st_crs               
#> [53] st_crs<-              st_difference        
#> [55] st_drop_geometry      st_filter            
#> [57] st_geometry           st_geometry<-        
#> [59] st_inscribed_circle   st_interpolate_aw    
#> [61] st_intersection       st_intersects        
#> [63] st_is_valid           st_is                
#> [65] st_join               st_line_merge        
#> [67] st_m_range            st_make_valid        
#> [69] st_nearest_points     st_node              
#> [71] st_normalize          st_point_on_surface  
#> [73] st_polygonize         st_precision         
#> [75] st_reverse            st_sample            
#> [77] st_segmentize         st_set_precision     
#> [79] st_shift_longitude    st_simplify          
#> [81] st_snap               st_sym_difference    
#> [83] st_transform          st_triangulate       
#> [85] st_union              st_voronoi           
#> [87] st_wrap_dateline      st_write             
#> [89] st_z_range            st_zm                
#> [91] summarise             transform            
#> [93] transmute             ungroup              
#> see '?methods' for accessing help and source code

Seguidamente, se describen y ejemplifican algunos de los métodos básicos de la clase sf.

13.5.3.2.1 st_read() - lectura de datos

El método st_read() lee datos vectoriales de una fuente en formato geoespacial (ej. shapefiles, archivos GeoJSON, bases de datos geoespaciales) y los recupera en un objeto sf.

En el siguiente bloque de código en R, se utiliza el método st_read() para leer un archivo GeoJSON con los polígonos de las provincias de Costa Rica. Este archivo proviene de un geoservicio de tipo Web Feature Service (WFS) publicado por el Instituto Geográfico Nacional (IGN).

# Lectura de una capa vectorial (GeoJSON) de provincias de Costa Rica
provincias <-
  st_read(
    "datos/ign/delimitacion-territorial-administrativa/provincias.geojson",
    quiet = TRUE # para evitar el despliegue de mensajes
  )

st_read() también puede crear objetos sf a partir de archivos de texto. Esta variante se utiliza principalmente cuando el archivo contiene coordenadas correspondientes a geometrías de puntos.

En el siguiente bloque de código, se utiliza st_read() para leer un archivo CSV con registros de presencia de felinos de Costa Rica, el cual contiene dos columnas llamadas decimalLongitude y decimalLatitude correspondientes a la longitud decimal y latitud decimal en las que fue observado cada felino. Este archivo proviene de una consulta al portal de datos de la Infraestructura Mundial de Información en Biodiversidad (GBIF).

# Lectura de un archivo CSV con registros de presencia de felinos en Costa Rica
felinos <-
  st_read(
    "datos/gbif/felinos.csv",
    options = c(
      "X_POSSIBLE_NAMES=decimalLongitude", # columna de longitud decimal
      "Y_POSSIBLE_NAMES=decimalLatitude"   # columna de latitud decimal
    ),
    quiet = TRUE
  )

Tanto provincias como felinos son objetos de la clase sf (y además de data.frame).

# Clase del objeto provincias
class(provincias)
#> [1] "sf"         "data.frame"

# Clase del objeto felinos
class(felinos)
#> [1] "sf"         "data.frame"

Al escribirse el nombre de un objeto sf en la consola de R, se despliega información general sobre este.

# Información general sobre el objeto provincias
provincias
#> Simple feature collection with 7 features and 7 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 156152 ymin: 608833.8 xmax: 658879.5 ymax: 1241118
#> Projected CRS: CR05 / CRTM05
#>                  gml_id cod_catalo cod_provin  provincia
#> 1 limiteprovincial_5k.1     160103          6 Puntarenas
#> 2 limiteprovincial_5k.2     160103          1   San José
#> 3 limiteprovincial_5k.3     160103          7      Limón
#> 4 limiteprovincial_5k.4     160103          3    Cartago
#> 5 limiteprovincial_5k.5     160103          2   Alajuela
#> 6 limiteprovincial_5k.6     160103          5 Guanacaste
#> 7 limiteprovincial_5k.7     160103          4    Heredia
#>                                                                                                                                                                                                                  ori_toponi
#> 1 En documento de 1720, se menciona la llegada del pirata Chipperton a la zona, en el cual aparece la descripcíon referente a una embarcación pequeña en la Punta de Arena, adoptando con el tiempo el nombre de Puntarenas
#> 2                                                                                                                                              Se remonta a la creación de la ermita dedicada al Patriarca San José en 1737
#> 3                                                                                         El origen del nombre de la provincia se remonta a 1852, cuando por primera vez se cita en un documento oficial el puerto de Limón
#> 4                                         Don Juan Vázques de Coronado escogió el sitio en el valle del Guarco para trasladar a la ciudad de Garcimuños, en 1563, bautizando al nuevo asentamiento con el nombre de Cartago
#> 5                                                                                                                  Se remonta al paraje llamado La Lajuela que por primera vez se cita en los Protocolos de Cartago de 1657
#> 6        En alegoria a un frondoso árbol de Guanacaste ubicado en la intersección de los caminos que se dirigían a Nicoya, Bagaces y Rivas, en lo que hoy día es el parque de Liberia. Esta referencia data del siglo XVIII
#> 7                                        En correspondiencia al Presidente  de la Real Audiencia de Guatemala, Capitán General don Alonso Fernández de Heredia, de la Inmaculada Concepción de Cubujuquí a Villa de Heredia
#>       area  version                       geometry
#> 1 11298.51 20201222 MULTIPOLYGON (((159917.7 60...
#> 2  4969.73 20201222 MULTIPOLYGON (((507877.9 11...
#> 3  9176.96 20201222 MULTIPOLYGON (((535996.1 11...
#> 4  3093.23 20201222 MULTIPOLYGON (((520206.6 11...
#> 5  9772.27 20201222 MULTIPOLYGON (((473981.3 11...
#> 6 10196.32 20201222 MULTIPOLYGON (((380807.3 11...
#> 7  2663.46 20201222 MULTIPOLYGON (((482356.4 11...
13.5.3.2.2 st_crs() y st_transform() - manejo de sistemas de coordenadas

El método st_crs() retorna el CRS de un objeto sf.

# Despliegue del CRS del objeto provincias
st_crs(provincias)
#> Coordinate Reference System:
#>   User input: CR05 / CRTM05 
#>   wkt:
#> PROJCRS["CR05 / CRTM05",
#>     BASEGEOGCRS["CR05",
#>         DATUM["Costa Rica 2005",
#>             ELLIPSOID["WGS 84",6378137,298.257223563,
#>                 LENGTHUNIT["metre",1]]],
#>         PRIMEM["Greenwich",0,
#>             ANGLEUNIT["degree",0.0174532925199433]],
#>         ID["EPSG",5365]],
#>     CONVERSION["Costa Rica TM 2005",
#>         METHOD["Transverse Mercator",
#>             ID["EPSG",9807]],
#>         PARAMETER["Latitude of natural origin",0,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8801]],
#>         PARAMETER["Longitude of natural origin",-84,
#>             ANGLEUNIT["degree",0.0174532925199433],
#>             ID["EPSG",8802]],
#>         PARAMETER["Scale factor at natural origin",0.9999,
#>             SCALEUNIT["unity",1],
#>             ID["EPSG",8805]],
#>         PARAMETER["False easting",500000,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8806]],
#>         PARAMETER["False northing",0,
#>             LENGTHUNIT["metre",1],
#>             ID["EPSG",8807]]],
#>     CS[Cartesian,2],
#>         AXIS["northing (N)",north,
#>             ORDER[1],
#>             LENGTHUNIT["metre",1]],
#>         AXIS["easting (E)",east,
#>             ORDER[2],
#>             LENGTHUNIT["metre",1]],
#>     USAGE[
#>         SCOPE["Cadastre, engineering survey, topographic mapping (large and medium scale)."],
#>         AREA["Costa Rica - onshore and offshore east of 86°30'W."],
#>         BBOX[2.21,-86.5,11.77,-81.43]],
#>     ID["EPSG",5367]]

# Despliegue CRS del objeto felinos
st_crs(felinos)
#> Coordinate Reference System: NA

st_crs() también puede asignar un CRS a un objeto sf que no lo tiene.

# Asignación de un CRS al objeto felinos
st_crs(felinos) <- 4326

El método st_transform() transforma un objeto sf a un nuevo CRS.

# Transformación del CRS del objeto provincias a WGS84 (EPSG = 4326)
provincias <-
  provincias %>%
  st_transform(4326)
13.5.3.2.3 plot() - mapeo

El método plot() grafica objetos sf en un mapa.

# Mapeo de las geometrías del objeto provincias
plot(provincias$geometry)

# Mapeo con argumentos adicionales de plot()
plot(
  provincias$geometry,
  extent = st_bbox(c(xmin = -86.0, xmax = -82.3, ymin = 8.0, ymax = 11.3)),
  main = "Provincias de Costa Rica",
  axes = TRUE,
  graticule = TRUE
)

Los argumentos reset y add de plot() permiten generar un mapa con varias capas.

# Primera capa del mapa
plot(
  provincias$geometry,
  extent = st_bbox(c(xmin = -86.0, xmax = -82.3, ymin = 8.0, ymax = 11.3)),
  main = "Registros de presencia de felinos en Costa Rica",
  axes = TRUE,
  graticule = TRUE,
  reset = FALSE
)

# Segunda capa
plot(felinos$geometry,
     add = TRUE,     
     pch = 16,
     col = "blue")

Para conocer los valores del argumento pch, puede consultar R plot pch symbols.

13.5.3.2.4 st_write() - escritura de datos

El método st_write() guarda en el disco un objeto sf en los diferentes formatos vectoriales de GDAL.

# Especificación del directorio de trabajo (debe utilizarse una ruta existente)
setwd("/home/mfvargas")

# Escritura del objeto provincias en formato GPKG
provincias %>%
  st_write("provincias.gpkg")

# Escritura del objeto felinos en formato KML
felinos %>%
  st_write("felinos.kml")
13.5.3.2.5 Otros
# Cantidad de filas de un objeto sf
nrow(provincias)
#> [1] 7

# Cantidad de columnas de un objeto sf
ncol(provincias)
#> [1] 8

# Resumen de la columna de geometría
summary(provincias$geometry)
#>  MULTIPOLYGON     epsg:4326 +proj=long... 
#>             7             0             0

13.5.4 Mapeo de objetos sf con otros paquetes

13.5.4.1 leaflet

El paquete leaflet genera mapas interactivos en lenguaje de marcado de hipertexto (HTML), el lenguaje de marcado utilizado para desarrollar páginas web.

13.5.4.1.1 Instalación y carga
# Instalación de leaflet
install.packages("leaflet")
# Carga de leaflet
library(leaflet)
13.5.4.1.2 Ejemplos

Antes de mostrar la capa de provincias en leaflet, se simplifica para evitar un error que se produce al incluir el mapa en un documento bookdown. Es posible que este paso sea innecesario en otros casos, por lo que podría omitirse.

provincias <-
  provincias %>%
  st_transform(5367) %>%
  st_simplify(dTolerance = 25) %>%
  st_transform(4326)

El siguiente bloque de código muestra las capas de provincias (polígonos) y de registros de presencia de felinos (puntos) sobre una capa base de OpenStreetMap (OSM).

# Mapa leaflet básico de provincias y registros de presencia de felinos
leaflet() %>%
  addTiles() %>% # capa base de OSM
  addPolygons( # capa de provincias (polígonos)
    data = provincias,
    color = "black",
    fillColor = "transparent",
    stroke = TRUE,
    weight = 1.0,
  ) %>%  
  addCircleMarkers( # capa de registros de presencia (puntos)
    data = felinos,
    stroke = F,
    radius = 4,
    fillColor = 'blue',
    fillOpacity = 1
  )

13.6 Datos raster

13.6.1 El modelo 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 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 13.2.

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)

Figure 13.2: El modelo raster: (A) ID de las celdas, (B) valores de las celdas, (C) mapa raster de colores. Imagen de Robin Lovelace et al.

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 13.3.

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

Figure 13.3: Ejemplos de mapas raster continuos y categóricos. Imagen de Robin Lovelace et al.

13.6.2 El paquete terra

El paquete terra implementa un conjunto de funciones para la lectura, escritura, manipulación, análisis y modelado de datos raster y vectoriales. Implementa la clase SpatRaster para manejar los objetos raster.

13.6.2.1 Instalación y carga

# Instalación de terra
install.packages("terra")
# Carga de terra
library(terra)

13.6.2.2 Métodos

La función help() presenta la documentación del paquete terra, incluyendo sus métodos.

# Ayuda sobre el paquete terra
help("terra-package")

Seguidamente, se describen y ejemplifican algunos de los métodos básicos para manejo de datos raster del paquete terra.

13.6.2.2.1 rast() - lectura de datos

El método rast() lee datos raster.

En el siguiente bloque de código en R, se utiliza el método rast()` para leer un archivo GeoTIFF correspondiente a la altitud de Costa Rica. Este archivo proviene de WorldClim, un conjunto de capas climáticas disponibles en varias resoluciones espaciales.

# Lectura de una capa raster de altitud
altitud <-
  rast(
    "datos/worldclim/altitud.tif"
  )

altitud es un objeto de la clase SpatRaster.

# Clase del objeto altitud
class(altitud)
#> [1] "SpatRaster"
#> attr(,"package")
#> [1] "terra"

Al escribirse el nombre de un objeto SpatRaster en la consola de R, se despliega información general sobre este.

# Información general sobre el objeto altitud
altitud
#> class       : SpatRaster 
#> dimensions  : 687, 546, 1  (nrow, ncol, nlyr)
#> resolution  : 0.008333333, 0.008333333  (x, y)
#> extent      : -87.10189, -82.55189, 5.494651, 11.21965  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : altitud.tif 
#> name        : altitud
13.6.2.2.2 crs() y project() - manejo de sistemas de coordenadas

El método crs() retorna el CRS de un objeto SpatRaster.

# CRS del objeto altitud
crs(altitud)
#> [1] "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"geodetic latitude (Lat)\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"geodetic longitude (Lon)\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    ID[\"EPSG\",4326]]"

crs() también puede asignar un CRS a un objeto SpatRaster.

# Asignación de un CRS a una copia del objeto altitud
altitud_crtm05 <- altitud
crs(altitud_crtm05) <- "EPSG:5367"

# Consulta
crs(altitud_crtm05)
#> [1] "PROJCRS[\"CR05 / CRTM05\",\n    BASEGEOGCRS[\"CR05\",\n        DATUM[\"Costa Rica 2005\",\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",5365]],\n    CONVERSION[\"Costa Rica TM 2005\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-84,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9999,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"northing (N)\",north,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"easting (E)\",east,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Cadastre, engineering survey, topographic mapping (large and medium scale).\"],\n        AREA[\"Costa Rica - onshore and offshore east of 86°30'W.\"],\n        BBOX[2.21,-86.5,11.77,-81.43]],\n    ID[\"EPSG\",5367]]"

El método project() reproyecta un objeto SpatRaster a un nuevo CRS.

# Transformación del CRS del objeto altitud
altitud_utm17N <-
  altitud %>%
  project("EPSG:8910")

# Consulta
crs(altitud_utm17N)
#> [1] "PROJCRS[\"CR-SIRGAS / UTM zone 17N\",\n    BASEGEOGCRS[\"CR-SIRGAS\",\n        DATUM[\"CR-SIRGAS\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",8907]],\n    CONVERSION[\"UTM zone 17N\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-81,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Oil and gas exploration and production.\"],\n        AREA[\"Costa Rica - offshore - Caribbean sea.\"],\n        BBOX[9.6,-83.6,11.77,-81.43]],\n    ID[\"EPSG\",8910]]"
13.6.2.2.3 plot() - mapeo

El método plot() muestra objetos SpatRaster en un mapa.

# Mapa de la capa de altitud
plot(altitud)

Los argumentos reset y add de plot() permiten generar un mapa con varias capas.

# Primera capa del mapa (raster)
plot(
  altitud,
  main = "Registros de presencia de felinos en Costa Rica",
  axes = TRUE,
  reset = FALSE
)

# Segunda capa (vectorial)
plot(felinos$geometry,
     add = TRUE,     
     pch = 16,
     col = "blue")
13.6.2.2.4 writeRaster() - escritura de datos

La función writeRaster() guarda en el disco un objeto SpatRaster en los diferentes formatos raster de GDAL.

# Especificación del directorio de trabajo (debe utilizarse una ruta existente)
setwd("/home/mfvargas")

# Escritura del objeto altitud
altitud %>%
  writeRaster("altitud.asc")
13.6.2.2.5 Otros
# Cantidad de filas de un objeto SpatRaster
nrow(altitud)
#> [1] 687

# Cantidad de columnas de un objeto SpatRaster
ncol(altitud)
#> [1] 546

# Resolución de un objeto SpatRaster
res(altitud)
#> [1] 0.008333333 0.008333333

13.6.3 Mapeo de objetos SpatRaster con otros paquetes

13.6.3.1 leaflet

A la fecha de escritura de este documento (2022-06-05), el paquete leaflet no cuenta con soporte para objetos de tipo SpatRaster, pero sí para objetos de tipo RasterLayer del paquete raster (antecesor de terra). Por este motivo, se utiliza el método raster() de raster para realizar la conversión de SpatRaster a RasterLayer y así mapear el objeto altitud con el método addRasterImage() de leaflet. Para más información sobre el tema de conversiones entre las diferentes clases espaciales de R, puede consultarse Conversions between different spatial classes in R.

Entonces, primero se instala y carga el paquete raster. También debe cargarse e instalarse el paquete rgdal, como complemento de raster.

# Instalación de raster y rgdal
install.packages("raster")
install.packages("rgdal")
# Carga de raster y rgdal
library(raster)
library(rgdal)

Se muestran en un mapa leaflet las capas de altitud (raster), provincias (polígonos) y registros de presencia de felinos (puntos) sobre una capa base de OpenStreetMap (OSM).

# Mapa leaflet básico con capas de altitud, provincias y registros de presencia de felinos
leaflet() %>%
  addTiles() %>%
  addRasterImage( # capa raster
    raster(altitud), # conversión de SpatRaster a RasterLayer 
    opacity = 0.6
  ) %>%    
  addPolygons(
    data = provincias,
    color = "black",
    fillColor = "transparent",
    stroke = TRUE,
    weight = 1.0,
  ) %>%
  addCircleMarkers(
    data = felinos,
    stroke = F,
    radius = 4,
    fillColor = 'blue',
    fillOpacity = 1
  )

13.7 Ejercicios

En un documento R Markdown incluya:

  1. Una tabla DT con las columnas del conjunto de datos de felinos de Costa Rica correspondientes a especie (species), fecha (eventDate), provincia (stateProvince), localidad (locality), longitud (decimalLongitude) y latitud (decimalLatitude).
    • Las columnas deben, si es necesario, convertirse al tipo de datos adecuado (ej. Date).
    • Los encabezados de las columnas en la tabla deben desplegarse en español, pero no deben alterarse los nombres de las columnas. Sugerencia: utilice el argumento colnames de la función datatable().
    • Los controles de la tabla deben estar en español.
  2. Un gráfico de barras, generado con gglot2 y traducido a plotly con ggplotly(), que muestre la cantidad de registros para cada especie de felinos.
    • Las barras deben estar ordenadas de mayor a menor.
    • Todos los controles y etiquetas del gráfico deben estar en español.
  3. Un mapa leaflet con las siguientes capas:
    • Capa base de OpenStreetMap (OSM).
    • Provincias de Costa Rica (polígonos).
    • Registros de presencia de felinos (puntos).

Publique el documento como un sitio web en GitHub Pages.

13.8 Recursos de interés

Bivand, R. (2022). CRAN Task View: Analysis of Spatial Data. https://CRAN.R-project.org/view=Spatial

R-Ladies Madrid. (2021). R-Ladies Madrid (español)—Analiza datos espaciales—Stephanie Orellana. https://www.youtube.com/watch?v=59tO2ARvVVU