Capítulo 6 Datos espaciales

6.1 Introducción

En este capitulo se presenta una breve introducción al uso de R como ambiente para trabajar datos espaciales. Los paquetes a utilizar son:

library(sp)
library(sf)
library(ggspatial)
library(raster)
library(terra)
library(stars)
library(viridis)
library(rgeos)
library(rgdal)
library(mapview)
library(RColorBrewer)
library(ggrepel)
library(rio)
library(tmap)
library(tidymodels)
library(tidyverse)
library(rayshader)

6.2 Paquetes para datos espaciales

La comunidad de R ha desarrollado una gran variedad de paquetes para datos espaciales, algunos de los cuales han ido evolucionando y facilitando la manipulación y presentación de los mismos.

Dentro de los paquetes más usados están:

sf (Pebesma, 2020 a): Para datos vectoriales dentro de la filosofía tidyverse
sp (Pebesma & Bivand, 2020 b): El predecesor de sf para datos vectoriales y en grilla (no exactamente raster)
terra y raster (Hijmans, 2020 b, 2020 a): Para datos raster
stars (Pebesma, 2020 b): El candidato a suceder a raster
rgeos y rgdal (Bivand et al., 2019; Bivand & Rundel, 2019): Interfaces para librerías básicas de manipulación de datos espaciales
ggplot2 y tmap (Tennekes, 2019): Creación de mapas estáticos (tmap puede crear mapas interactivos)
mapview y leaflet (Appelhans et al., 2020; Cheng et al., 2018): Creación de mapas interactivos

Otro montón de paquetes brindan funciones adicionales y funciones para tareas especificas.

6.3 Sistemas de Referencias de Coordenadas (CRS)

Los datos espaciales tienen por lo general un sistema de coordenadas asociado (geográficas, lambert, UTM, etc.). Estos sistemas se pueden identificar por medio de códigos EPSG que han sido estandarizados para uso general.

Como ejemplo, las coordenadas geográficas en grados (WGS84) tiene el código 4326, las coordenadas geográficas en metros (WGS84/pseudo-mercator; son las utilizadas por GoogleMaps y otros) tiene el codigo 3857. Las coordenadas para Costa Rica con la proyección CRTM05 tienen el codigo 5367, mientras que para la proyección Lambert Norte tienen el codigo 5456. Todos estos se pueden obtener por medio del sitio web epsg.io o en QGIS a la hora de buscar sistemas de coordenadas ahí aparece el codigo.

6.4 Importar datos

6.4.1 Desde archivos de texto

Datos en archivos de texto (.txt, .csv, etc. y principalmente para datos puntuales) se pueden importar de manera convencional con rio::import.

datos <- rio::import("data/BroomsBarn.txt", setclass = 'tibble') %>% 
  mutate(x = x*40, y = y*40, logK = log(K), logP = log(P))

Posteriormente pueden convertirse a datos espaciales (st_as_sf), indicando la posición o el nombre de las columnas con las coordenadas X,Y y asignándole un CRS.

datos_sf = st_as_sf(datos, coords = 1:2, crs = NA, remove = F)

La ventaja de sf es que permite manipular los dato asociados al objeto espacial haciendo uso de los verbos del tidyverse.

datos_sf %>% 
  filter(pH > 8)

## # A tibble: 163 x 8
##        x     y     K  logK    pH     P  logP   geometry
##    <dbl> <dbl> <int> <dbl> <dbl> <dbl> <dbl>    <POINT>
##  1    80  1240    14  2.64   8.2   6.6  1.89  (80 1240)
##  2   120    40    28  3.33   8.2   4.2  1.44   (120 40)
##  3   120    80    26  3.26   8.2   7.5  2.01   (120 80)
##  4   120   120    23  3.14   8.2   5.1  1.63  (120 120)
##  5   120   160    21  3.04   8.2   4    1.39  (120 160)
##  6   120   240    22  3.09   8.3   4.8  1.57  (120 240)
##  7   120   280    24  3.18   8.3   4.8  1.57  (120 280)
##  8   120   320    41  3.71   8.2   4    1.39  (120 320)
##  9   120  1240    15  2.71   8.2   8.4  2.13 (120 1240)
## 10   160    80    24  3.18   8.2   6.7  1.90   (160 80)
## # … with 153 more rows

Para algunas funciones todavía es necesario usar objetos sp por lo que se pueden crear estos a partir de los objetos sf.

datos_sp = as(datos_sf, 'Spatial')
coordnames(datos_sp) = c('X','Y')

6.4.2 Desde archivos espaciales

6.4.2.1 Shapefiles

De igual manera se pueden importar/leer directamente datos en muchos formatos espaciales (ejemplo: shapefiles, geopackage, raster, etc.). Las funciones para leer datos son: st_read y read_sf; la primera importa los datos como DataFrame, la segunda como tibble.

fallas = read_sf('data/fallas.shp')
geomorfo = read_sf('data/geomorfo.shp')

En este caso no reconocen los metadatos de la proyeccion por lo que se les puede asignar por medio de st_set_crs, donde el argumento necesario es el codigo EPSG.

fallas_ln = fallas %>% 
  st_set_crs(5456)
geomorfo_ln = geomorfo %>% 
  st_set_crs(5456)

Si se desean transformar a otro sistema se usa st_transform, donde, de nuevo, el argumento necesario es el codigo EPSG del sistema destino (En el caso del ejemplo a coordenadas geográficas).

fallas_geog = fallas_ln %>% 
  st_transform(4326)
geomorfo_geog = geomorfo_ln %>% 
  st_transform(4326)

6.4.2.2 Geopackage

Si se tiene un archivo .gpkg es necesario explorar primero las capas disponibles para luego importar los datos deseados. Esto se realiza con st_layers y la dirección del archivo geopackage.

st_layers(dsn = 'data/espaciales.gpkg')

## Driver: GPKG 
## Available layers:
##   layer_name geometry_type features fields
## 1     fallas   Line String     2028      7
## 2   geomorfo       Polygon      467      6

Una vez se han identificado las capas se pueden importar con cualquiera de los métodos deseados, donde el argumento a especificar es layer y el nombre de la capa a importar.

fallas2 = read_sf('data/espaciales.gpkg', layer = 'fallas')
geomorfo2 = read_sf('data/espaciales.gpkg', layer = 'geomorfo')

6.4.2.3 Raster

Para leer rasters de una banda se usa la función raster, para rasters multi-banda se puede usar brick. Para objetos stars es indiferente y se pueden leer por medio de read_stars. Para convertir de raster a stars se usa st_as_stars.

La lectura de un raster de una banda se ejemplifica con un modelo de elevación digital del Pacifico.

pacifico = raster('data/pacifico.tif')
pacifico

## class      : RasterLayer 
## dimensions : 300, 300, 90000  (nrow, ncol, ncell)
## resolution : 0.01661111, 0.01661111  (x, y)
## extent     : -87.99167, -83.00833, 5.008333, 9.991667  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : /Users/maximiliano/Documents/UCR/Docencia/Extras/R/bookdown/geolonum/data/pacifico.tif 
## names      : pacifico 
## values     : -4770, 3601  (min, max)

pacifico_stars = read_stars('data/pacifico.tif')
pacifico_stars = st_as_stars(pacifico)
pacifico_stars

## stars object with 2 dimensions and 1 attribute
## attribute(s):
##  pacifico.tif   
##  Min.   :-4770  
##  1st Qu.:-3144  
##  Median :-2417  
##  Mean   :-2155  
##  3rd Qu.:-1627  
##  Max.   : 3601  
## dimension(s):
##   from  to   offset      delta                       refsys point values    
## x    1 300 -87.9917  0.0166111 +proj=longlat +datum=WGS8... FALSE   NULL [x]
## y    1 300  9.99167 -0.0166111 +proj=longlat +datum=WGS8... FALSE   NULL [y]

La lectura de un raster multi-banda se ejemplifica con un geotiff que viene con el paquete stars.

sat_ras = brick(system.file('tif/L7_ETMs.tif', package = 'stars'))
sat_ras

## class      : RasterBrick 
## dimensions : 352, 349, 122848, 6  (nrow, ncol, ncell, nlayers)
## resolution : 28.5, 28.5  (x, y)
## extent     : 288776.3, 298722.8, 9110729, 9120761  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=25 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## source     : /Library/Frameworks/R.framework/Versions/3.6/Resources/library/stars/tif/L7_ETMs.tif 
## names      : L7_ETMs.1, L7_ETMs.2, L7_ETMs.3, L7_ETMs.4, L7_ETMs.5, L7_ETMs.6 
## min values :         0,         0,         0,         0,         0,         0 
## max values :       255,       255,       255,       255,       255,       255

sat_df = as.data.frame(sat_ras,xy=T)
sat_stars = read_stars(system.file('tif/L7_ETMs.tif', package = 'stars'))
sat_stars

## stars object with 3 dimensions and 1 attribute
## attribute(s):
##   L7_ETMs.tif    
##  Min.   :  1.00  
##  1st Qu.: 54.00  
##  Median : 69.00  
##  Mean   : 68.91  
##  3rd Qu.: 86.00  
##  Max.   :255.00  
## dimension(s):
##      from  to  offset delta                       refsys point values    
## x       1 349  288776  28.5 PROJCS["UTM Zone 25, Sout... FALSE   NULL [x]
## y       1 352 9120761 -28.5 PROJCS["UTM Zone 25, Sout... FALSE   NULL [y]
## band    1   6      NA    NA                           NA    NA   NULL

6.5 Exportar datos

6.5.1 Vectoriales

Para exportar datos vectoriales se puede usar st_write, donde se define el objeto espacial a exportar y el tipo de archivo a generar (.shp, .gpkg, etc.). Las opciones layer_options = 'OVERWRITE=YES' y delete_layer = T permiten reescribir un objeto si ya se encontraba presente en la carpeta destino.

st_write(fallas_geog, dsn = 'data/espaciales_geog.gpkg', 
         layer = 'fallas', layer_options = 'OVERWRITE=YES', 
         quiet = T, delete_layer = T)
st_write(geomorfo_geog, dsn = 'data/espaciales_geog.gpkg', 
         layer = 'geomorfo', layer_options = 'OVERWRITE=YES', 
         quiet = T, delete_layer = T)

Se puede revisar que el objeto haya sido creado apropiadamente, de nuevo usando st_layers.

st_layers(dsn = 'data/espaciales_geog.gpkg')

## Driver: GPKG 
## Available layers:
##   layer_name geometry_type features fields
## 1     fallas   Line String     2028      7
## 2   geomorfo       Polygon      467      6

6.5.2 Raster

Dependiendo de si el objeto es raster o stars, se deberá usar la función respectiva, writeRaster o write_stars. en el caso de writeRaster se tiene que especificar el formato (para esto se puede consultar la ayuda de la función).

writeRaster(sat_ras, 'data/imagen_satelite.tif', 
            format="GTiff", overwrite=TRUE)
write_stars(sat_stars, 'data/imagen_satelite_stars.tif')

6.6 Mapas

Existen diferentes formas de graficar datos espaciales, aquí se presentan las más usadas, donde se recomienda ggplot2 o tmap. El primero ya que se esta familiarizado con el funcionamiento del mismo y es nada más agregar un par de funciones especificas para datos espaciales; el segundo ya que es especifico para datos espaciales siguiendo una idea similar a ggplot2.

6.6.1 Estáticos

6.6.1.1 básico

Todos los objetos espaciales tienen un método de ploteo básico (plot). Si el objeto tiene más de 1 atributo va a plotear todos o los que pueda. Para evitar esto se puede especificar cual atributo se quiere plotear. En general esto se puede usar más para una visualización rápida pero no para mapas finales. Se muestran diferentes opciones para visualizar los diferentes tipos de datos, pero no se muestra el mapa final por ser muy básico.

plot(fallas_ln)

plot(datos_sf[,'K'])

plot(geomorfo_ln[,'CODIGO'])

plot(pacifico)

plot(pacifico_stars)

6.6.1.2 ggplot

El paquete ggplot tiene geometrías especificas para objetos sf (vector) y stars (raster o grillas), por lo que facilita la creación de mapas en un ambiente ya conocido.

La forma más básica de crear un mapa vectorial usando ggplot es usando geom_sf y especificando el argumento data que corresponde con el objeto espacial, pero por defecto despliega coordenadas geográficas (Figura 6.1).

ggplot() + 
  geom_sf(data = fallas_ln)

Figura 6.1: Mapa básico en ggplot2

Para que el mapa se despliegue en las coordenadas del objeto y no geográficas, hay que cambiar el argumento datum en coord_sf al CRS del archivo deseado (Figura 6.2).

ggplot() + 
  geom_sf(data = fallas_ln) + 
  coord_sf(datum = st_crs(fallas_ln))

Figura 6.2: Mapa en ggplot2 con el datum modificado

El paquete ggspatial (Dunnington, 2018) ofrece algunas capas adicionales especificas (elementos cartográficos) para datos espaciales, con las cuales se pueden agregar una escala (annotation_scale) y el norte (annotation_north_arrow), ademas de poder agregar un mapa de fondo (annotation_map_tile), donde zoom define el nivel de detalle, a menor zoom menor detalle (Figura 6.3).

ggplot() + 
  annotation_map_tile(zoom = 8) +
  geom_sf(data = fallas_ln) + 
  coord_sf(datum = st_crs(fallas_ln)) + 
  annotation_scale(location = 'bl') + 
  annotation_north_arrow(location = 'tr',
                         height = unit(.75, "cm"), 
                         width = unit(.75, "cm"))

Figura 6.3: Mapa en ggplot con elementos cartográficos (escala y norte)

Para datos puntuales (Figura 6.4) o de polígonos (Figura 6.5) se puede especificar una columna de los datos, por la cual colorear los puntos.

ggplot() + 
  geom_sf(data = datos_sf, aes(col = K), size = 3, alpha = 0.6) + 
  scale_color_viridis_c()

Figura 6.4: Mapa de puntos en ggplot coloreados por la variable “K”

ggplot() + 
  geom_sf(aes(fill = as_factor(FORMA)), data = geomorfo_ln) + 
  coord_sf(datum = st_crs(geomorfo_ln)) + 
  scale_fill_brewer(palette = 'Set3')

Figura 6.5: Mapa de polígonos en ggplot coloreados por la variable “FORMA”

La forma de crear un mapa a partir de un objeto stars (grilla de una banda) es por medio de geom_stars, donde, de nuevo, se especifica el objeto espacial en el argumento data. En este caso automáticamente se aplica al relleno (fill) la variable del objeto. Ademas, es necesario agregar coord_equal para que los ejes de coordenadas tengan las misma escala (Figura 6.6).

ggplot() + 
  geom_stars(data = pacifico_stars) + 
  scale_fill_gradient(low = 'black', high = 'white', na.value = 'white') +
  coord_equal()

Figura 6.6: Mapa de objeto stars en ggplot. Se especifica el relleno como un gradiente y se aplica la misma escala a ámbos ejes

Para imágenes multi-banda es necesario transforma el objeto espacial en una tabla (sat_df = as.data.frame(sat_ras,xy=T)) y graficar el relleno (fill) usando la función rgb, donde se especifica el nombre de las bandas que corresponde con cada uno de los canales rgb, y se debe definir maxColorValue=255 (Figura 6.7). Como lo que se va a graficar es una grilla se usa la función geom_tile, y de nuevo es necesario definir los ejes de coordenadas iguales con coord_equal.

ggplot(sat_df, aes(x, y, 
                   fill=rgb(red = L7_ETMs.3, 
                            green = L7_ETMs.2, 
                            blue = L7_ETMs.1, 
                            maxColorValue = 255))) + 
  geom_tile() + 
  scale_fill_identity() + 
  coord_equal() + 
  scale_x_continuous(expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0)) +
  theme(axis.text.y = element_text(angle = 90, hjust = .5))

Figura 6.7: Mapa de color verdadero para una imagen satelital multi-banda

6.6.1.3 tmap

El paquete tmap es una opción especifica para datos espaciales y sigue la misma ideología de ggplot al trabajar en capas. Por defecto no despliega la grilla de coordenadas, hay que agregárselas con tm_grid. Para revertir los colores en la paleta se le agrega un menos (-) antes del nombre, y para colocar la leyenda fuera del área del mapa se debe usar tm_layout(legend.outside = T). Existen opciones globales para los mapas generados con tmap, estas se pueden modificar con tmao_options, y uno de los argumentos es la leyenda afuera, así no hay que aplicarlo en cada mapa por separado.

tmap_options(legend.outside = T)

La estructura básica de cualquier gráfico es:

tm_shape(shp = <DATA>) +
  <tm_function>(col = <'VARIABLE'>, palette = '', 
                style = '', n = 5) +
  <tm_layout> +
  <tm_xlab> + 
  <tm_ylab> + 
  <tm_grid> +
  <tm_scale_bar> + 
  <tm_compass> +

Se pueden agregar diferentes objetos espaciales agregando diferentes tm_shape. A diferencia de ggplot y geom_sf que era una función para todo objeto vectorial, es necesario definir el tipo de objeto con las diferentes funciones tm_*. Lo que se va a hacer es recrear los mapas de la sección anterior pero con tmap.

Un mapa básico se observa en la Figura 6.8, donde se grafican las líneass de falla.

tm_shape(fallas_ln) + 
  tm_lines()

Figura 6.8: Mapa básico con tmap

El mapa de puntos se crea por medio de tm_dots. Por defecto tmap discretiza los datos (Figura 6.9), pero si se quiere representar como un gradiente (Figura 6.10) es necesario definir style = 'cont', para continuo.

tm_shape(datos_sf) +
  tm_dots('K', size = .3, palette = 'viridis')

Figura 6.9: Mapa de puntos con tmap coloreado de acuerdo a una variable

tm_shape(datos_sf) +
  tm_dots('K', size = .3, palette = 'viridis', style = 'cont') + 
  tm_layout(legend.outside = T)

Figura 6.10: Mapa de puntos con tmap coloreado de acuerdo a una variable, con la leyenda tipo gradiente

En la Figura 6.11 se observa un mapa de polígonos, donde ademas se agrega la grilla de coordenadas pero sin las líneas.

geomorfo_ln %>% 
  mutate(FORMA = as_factor(FORMA)) %>% 
tm_shape() + 
  tm_polygons('FORMA') + 
  tm_grid(lines = F)

Figura 6.11: Mapa de poligonos con tmap, agregando la grilla de coordenadas sin la líneas internas

Para graficar objetos raster de una banda se usa tm_raster (Figura 6.12), donde hay que definir la paleta de colores a usar. Si se desea invertir simplemente se le agrega un - en frente del nombre.

tm_shape(pacifico) + 
  tm_raster(palette = '-Greys', style = 'cont')

Figura 6.12: Mapa raster con tmap

Para ver las opciones de las paletas de color se puede usar:

tmaptools::palette_explorer()

Para graficar objetos raster multi-banda se debe usar tm_rgb (Figura 6.13) y especificar el orden de las bandas que corresponde con cada uno de los canales rgb.

tm_shape(sat_ras) + 
  tm_rgb(r = 3, g = 2, b = 1)

Figura 6.13: Mapa de color verdadero con tmap

Finalmente, tmap trae incorporadas funciones para agregar escala (tm_scale_bar) y norte (tm_compass) (Figura 6.14).

tm_shape(fallas_ln) + 
  tm_lines() + 
  tm_scale_bar(width = .2,position = c('left','bottom')) + 
  tm_compass(position = c('right','top'))

Figura 6.14: Mapa incorporando elementos de escala y norte

Para salvar un mapa de tmap como imagen se usa tmap_save. Idealmente hay que guardar el mapa en un objeto.

mapa1 = tm_shape(fallas_ln) + 
  tm_lines() + 
  tm_scale_bar(width = .2,position = c('left','bottom')) + 
  tm_compass(position = c('right','top'))

tmap_save(mapa1, 'figures/mapa_fallas.png', dpi = 300)

6.6.2 Dinámicos

Una de la opciones con mapas interactivos es que se pueden definir diferentes mapas de fondo. Dentro de los más usados están los que se muestran a continuación:

mybasemaps = c('CartoDB.Positron', 'OpenStreetMap', 'OpenTopoMap',
               'Esri.WorldImagery', 'Esri.WorldTopoMap', 'Esri.OceanBasemap')

De manera general se pueden asignar mapas de fondo para todos los mapas de una sesión definiéndolos en las opciones de tmap y mapview:

tmap_options(basemaps = mybasemaps)

mapviewOptions(basemaps = mybasemaps)

La Figura 6.15 muestra las diferentes opciones para manipular un mapa dinámico creado con mapview. Las mismas opciones y elementos se pueden agregar a un mapa de leaflet pero requieren ser especificadas explícitamente (haciendo más largo el codigo), mientras que mapview ya agrega muchos de estos por defecto.

Figura 6.15: Partes de un mapa dinámico creado con mapview

6.6.2.1 tmap

Este paquete ofrece la opción de visualizar un mapa estático como dinámico, cambiando el modo de visualización (tmap_mod(mode = c('plot','view'))), donde plot es para mapas estático y view para mapas dinámicos y traduce el mapa a un mapa leaflet.

Lo anterior cambia el modo para todos los mapas que se creen a partir de que esto se modifica. Para visualizar de manera dinámica un mapa guardado en un objeto se puede usar tmap_leaflet, que genera un mapa leaflet a partir de uno de tmap, pero solo para el mapa que se quiere (Figura 6.16). Para demostrar esto se genera el mapa dinámico a partir del mapa que se salvo anteriormente (mapa1), que corresponde con la Figura 6.8.

tmap_leaflet(mapa1)

Figura 6.16: Mapa interactivo leaflet a partir de un mapa estático tmap

6.6.2.2 mapview

Este paquete permite generar mapas interactivos de manera eficiente y rápida, pero no brinda la personalización de leaflet, el cual es mucho más complejo y para generar un mapa similar se requiere aproximadamente de 4 a 5 veces mas líneass de codigo.

La función básica es mapview donde se le pasa el objeto espacial. Por defecto le asigna un color único. El argumento layer_name es para definir el nombre que se quiere aparezca en la leyenda. Si se desean agregar diferentes capas simplemente se agregan funciones mapview con el operador + (Figura 6.17). mapview es inteligente en el sentido de que dependiendo del tipo de objeto (polígono, líneas, punto) lo va a graficar en el orden jerarquico visual (puntos por encima de líneass, líneass por encima de polígonos).

mapview(fallas_ln, layer.name = 'Fallas', color = 'red') +
  mapview(geomorfo_ln, layer.name = 'Geomorfologia')

Figura 6.17: Mapa con mapview, donde se pueden combinar objetos espaciales usando el operador +

Se pueden colorear objetos espaciales definiendo una columna de la tabla de atributos por medio del argumento zcol (Figura 6.18).

mapview(datos_sf, zcol = 'logK', layer.name = 'logK')

Figura 6.18: Mapa de puntos con mapview coloreados de acuerdo a una variable

Para graficar objetos raster de una banda simplemente se pasa el objeto a la función mapview (Figura 6.19). Para raster multi-banda se debe usar viewRGB y definir el orden de las bandas de acuerdo a los canales rgb (Figura 6.20).

mapview(pacifico)

Figura 6.19: Mapa raster con mapview

viewRGB(sat_ras, r = 3, g = 2, b = 1)

Figura 6.20: Mapa de color verdradero con mapview

Con mapview se pueden crear diferentes mapas y sincronizarlos por medio de sync del paquete leafsync (Simplemente se muestra el codigo).

m1 = mapview(franconia, zcol = 'district', 
             layer.name = 'Distrito')
m2 = mapview(breweries, legend = F)

leafsync::sync(m1, m2)

6.6.3 Modelos de sombras

Para generar modelos de sombras en 2D y 3D se usa el paquete rayshader (Morgan-Wall, 2019), el cual funciona a partir de objetos raster. Funciona con el pipe operator (%>%), por lo que es familiar a trabajar en el tidyverse.

Primero hay que pasar el raster a matriz.

pacifico_mat = raster_to_matrix(pacifico)

Para mejorar la apariencia del modelo de sombras es recomendado agregar diferentes sombras (‘shades’), las cuales son calculadas de diferente manera para resaltar diferentes aspectos del modelo.

zscale = 200
elmat = pacifico_mat
raymat = ray_shade(elmat, zscale = zscale)
ambmat = ambient_shade(elmat, zscale = zscale)
lambmat = lamb_shade(elmat, zscale = zscale)

Para generar un modelo de sombras en 2D, se empieza con la matriz de elevación, se le agregan las diferentes capas, y para mostrar el mapa se termina la cadena de comandos con plot_map (Figura 6.21).

elmat %>%
  sphere_shade(texture = "bw") %>%
  add_shadow(raymat, 0.5) %>%
  add_shadow(lambmat, 0.5) %>%
  add_shadow(ambmat, 0.5) %>%
  plot_map()

Figura 6.21: Modelo de sombras en 2D

Para generar el modelo en 3D (Figura 6.22), simplemente se termina la cadena de comandos con plot_3d. Esto abre una ventana interactiva donde se puede manipular el modelo, en este caso simplemente se toma una captura de la apariencia en 3D.

elmat %>%
  sphere_shade(texture = "imhof1") %>%
  add_shadow(raymat, 0.5) %>%
  add_shadow(lambmat, 0.5) %>%
  add_shadow(ambmat, 0.5) %>%
  plot_3d(elmat, zscale = zscale, water = T, 
          theta = 0, phi = 40, zoom = .5,
          windowsize = c(1000,600))
render_snapshot(clear=TRUE)

Figura 6.22: Modelo de sombras en 3D

Adicionalmente se puede crear una pequeña película o animación con render_movie (No se presenta aquí pero se muestra el codigo).

elmat %>%
  sphere_shade(texture = "imhof1") %>%
  add_shadow(raymat, 0.5) %>%
  add_shadow(lambmat, 0.5) %>%
  add_shadow(ambmat, 0.5) %>%
  plot_3d(elmat, zscale = zscale, water = T, 
          theta = 0, phi = 40, zoom = .5,
          windowsize = c(1000,600))
render_movie(filename = 'modelo3D')
rgl::rgl.close()

6.7 Recursos

Se presentan recursos a consultar para ahondar más en los temas presentados.

Geocomputatio with R Es un libro virtual donde se introducen conceptos de como usar R para análisis espacial y creación de mapas.

Spatial Data Sience

Intro to GIS and Spatial Analysis

RSpatial

Geoestadistica con R (El sitio web está en polaco, pero se puede traducir en Chrome)

Introducción a estadística con R (Capítulo 7)

Spatial workshop notes

Referencias

Appelhans, T., Detsch, F., Reudenbach, C., & Woellauer, S. (2020). mapview: Interactive Viewing of Spatial Data in R. https://CRAN.R-project.org/package=mapview

Bivand, R., Keitt, T., & Rowlingson, B. (2019). rgdal: Bindings for the ’Geospatial’ Data Abstraction Library. https://CRAN.R-project.org/package=rgdal

Bivand, R., & Rundel, C. (2019). rgeos: Interface to Geometry Engine - Open Source (’GEOS’). https://CRAN.R-project.org/package=rgeos

Cheng, J., Karambelkar, B., & Xie, Y. (2018). leaflet: Create Interactive Web Maps with the JavaScript ’Leaflet’ Library. https://CRAN.R-project.org/package=leaflet

Dunnington, D. (2018). ggspatial: Spatial Data Framework for ggplot2. https://CRAN.R-project.org/package=ggspatial

Hijmans, R. J. (2020a). raster: Geographic Data Analysis and Modeling. https://CRAN.R-project.org/package=raster

Hijmans, R. J. (2020b). terra: Classes and Methods for Spatial Data. https://CRAN.R-project.org/package=terra

Morgan-Wall, T. (2019). rayshader: Create and Visualize Hillshaded Maps from Elevation Matrices. https://CRAN.R-project.org/package=rayshader

Pebesma, E. (2020a). sf: Simple Features for R. https://CRAN.R-project.org/package=sf

Pebesma, E. (2020b). stars: Spatiotemporal Arrays, Raster and Vector Data Cubes. https://CRAN.R-project.org/package=stars

Pebesma, E., & Bivand, R. (2020b). sp: Classes and Methods for Spatial Data. https://CRAN.R-project.org/package=sp

Tennekes, M. (2019). tmap: Thematic Maps. https://CRAN.R-project.org/package=tmap