English-language version here
Hoy en día ya es posible y además relativamente fácil renderizar (anglicismo para describir el proceso de generar imágenes digitales a partir de un modelo tridimensional) mapas con perspectiva e iluminación usando R. Gracias a los paquetes rayshader y rayrender, desarrollados por Tyler Morgan-Wall y elevatr de Jeff Hollister, podemos combinar una imagen georeferenciada con un modelo de elevación digital (DEM, por sus siglas en inglés) para darle profundidad y ambiente a cualquier mapa que tengamos en formato digital.
El tuitero @flotsam siempre comparte ejemplos muy buenos, como éste que subió hace poco:
A 1972 NRCan (formerly the Department of Energy, Mines, and Resources) map of Montréal, Canada. Spent more time trying to decipher cryptic filenames to find what I want than actually rayshading it. :P#rayshader adventures, an #rstats tale pic.twitter.com/O756ddhs9G
— flotsam (@researchremora) August 2, 2021
El proceso para renderizar imagenes georeferencias está resumido en esta guía.
Hay muchos mapas ya georeferenciados en Internet, y además podemos usar programas de SIG para empatar cualquier mapa digital con su contexto espacial mediante puntos de control. Yo generalmente prefiero hacer mis propios mapa en ggplot2 y quería procesarlos en 3D con rayshader.
Cuando hacemos mapas en ggplot2 con los paquetes geoespaciales más populares (sf, stars, etc.), ya tenemos datos espaciales y no haría falta asignarle a mano coordenadas del mundo real a cada pixel del raster de nuestro mapa. Además, cada que trato de georeferenciar imágenes en qGIS pierdo la paciencia por tener que crear al menos seis puntos de control con el cursor y la mano toda temblorina.
En esta respuesta de StackOverflow, nos explican como extraer la información que describe los valores xy máximos y mínimos del panel de un objeto ggplot, para poder asignarle esta extensión a la versión raster que ya exportamos de nuestro objeto.
Para poder seguir este método, lo importante es asegurarse que el objeto de ggplot no tenga nada de bordes, para que los valores xy que definen su extensión empaten con los bordes de la imagen que exportamos con el fin de usarla en rayshader. Para lograrlo, usamos el (theme_nothing) que viene con el paquete cowplot de Claus Wilke. En resumen, ésto equivale a apagar casi todos los elementos de la figura y en no tener ninguna expansión de los ejes. Aquí podemos ver cuáles elementos se están modificando.
Probemos. Este ejemplo está basado en datos espaciales del portal de datos geográficos de la CONABIO y de la CONANP. Descargué los shapefiles de Áreas Naturales Protegidas y de Uso del suelo y vegetación (escala 1:250000, serie VI continuo nacional), y me enfoqué en la ANP Iztaccíhuatl-Popocatépetl.
Estos pasos cargan los archivos shp, recortan los datos a la ANP de interés, y agrupan (medio arbitrariamente) algunas de las categorías de uso de suelo y vegetación. Luego, podemos hacer un mapa simple pero colorido de los volcanes. Aquí el theme_nothing() le quita los márgenes y los rótulos de los ejes a la figura, y hay algunos argumentos de coord_sf() que son para que todo quede dentro del panel principal de la figura.
# paquetes
library(sf) # CRAN v0.9-8
library(dplyr) # CRAN v1.0.7
library(stringr) # CRAN v1.4.0
library(ggplot2) # CRAN v3.3.3
library(raster) # CRAN v3.4-13
library(stars) # [github::r-spatial/stars] v0.5-3
library(elevatr) # CRAN v0.4.1
library(rayshader) # [github::tylermorganwall/rayshader] v0.26.1
library(smoothr) # CRAN v0.2.2
# importar ANPs
anps <- st_read("YOUR-PATH-HERE/182ANP_Geo_ITRF08_Agosto_2020.shp") %>%
dplyr::select(NOMBRE)
# importar vegetación
veg <- st_read("YOUR-PATH-HERE/usv250s6gw.shp") %>%
st_transform(st_crs(anps))
# area de estudio
izta <- anps %>% filter(str_detect(NOMBRE, "Izt"))
# limites
Iztbbox <- st_buffer(izta, dist = 0.1) %>%
st_bbox() %>%
st_as_sfc()
# disolver vegetación
veggr <- veg %>%
group_by(CVE_UNION) %>%
summarize(.groups = "keep")
# recortar a los limites
Iztveg <- st_crop(veggr, Iztbbox)
# reclasificar vegetación
Iztvegrc <- Iztveg %>% mutate(vegclass = case_when(
str_detect(CVE_UNION, "^B") ~ "Bosque",
str_detect(CVE_UNION, "^VS") ~ "Vegetación secundaria",
str_detect(CVE_UNION, "^T") ~ "Agricultura (temporal)",
str_detect(CVE_UNION, "^R") ~ "Agricultura (Riego)",
str_detect(CVE_UNION, "DV") ~ "Sin vegetación",
str_detect(CVE_UNION, "^P") ~ "Pastizal",
str_detect(CVE_UNION, "VW") ~ "Pradera de alta montaña",
TRUE ~ CVE_UNION
))
# smooth
Iztvegrc <- smooth(Iztvegrc)
# exportando a tiff
ragg::agg_tiff(
filename = "izt.tiff",
res = 400, width = 3.88, height = 5.92, units = "in"
)
# figura
iztgg <-
ggplot() +
geom_sf(data = Iztvegrc, aes(fill = vegclass), size = 0) +
scale_fill_manual(values = c(
"#b5c99a", "#97a97c", "#8d99ae", "#26432F",
"#e8e8e4", "#cdeac0", "#B1D5BB",
"#cde5d7", "#adb6c4", "#9cc5a1",
"#ffc49b", "#eaeaea", "#f3e7e4"
), guide = FALSE) +
ggfx::with_outer_glow(geom_sf(data = izta, fill = "transparent", color = "white", size = 0.1),
expand = 3, colour = "white"
) +
coord_sf(expand = FALSE, clip = "on") +
cowplot::theme_nothing() +
theme(
panel.grid.major = element_line(color = "#14213d", size = 0.1, linetype = "dashed"),
panel.ontop = TRUE
)
iztgg
dev.off()Me tomó algunos intentos poder exportar el tiff sin que saliera espacio en blanco en los márgenes, pero cuando ésto ya queda podemos importar la imagen usando el paquete raster.
Revisemos ese tiff.
# creando un StackedRaster
stackedRaster <- raster::stack("izt.tiff")
La información geoespacial que necesitamos está contenida en el objeto ggplot, y la podemos extraer con ggplot_build. Uno de los componentes del objeto es una lista con los máximos y mínimos para x y para y, y éstos se los damos a la función extent para georeferenciar el mapa (también hay que definir la proyección). El raster georeferenciado se exporta con INT1U para que los valores queden entre 0 y 255, y PHOTOMETRIC=RGB para declarar cómo se interpretan los colores.
# extraer componentes Geoepaciales
lat_long <- ggplot_build(iztgg)$layout$panel_params[[1]][c("x_range", "y_range")]
# asignarle los valores al StackedRaster
raster::extent(stackedRaster) <- c(lat_long$x_range, lat_long$y_range)
raster::projection(stackedRaster) <- raster::crs(st_crs(Iztvegrc)$proj4string)
# exportar
writeRaster(stackedRaster, "iztGeoTiff.tif", options = "PHOTOMETRIC=RGB", datatype = "INT1U", overwrite = TRUE)Importamos de nuevo el geoTiff
grRast <- raster::brick("iztGeoTiff.tif")
grRastst <- raster::stack(grRast)Ahora ya podemos seguir las guías existentes para procesar imágenes con rayshader. Primero hay que transformar el raster en un array (arreglo de valores numéricos), y usar este objeto para descargar datos de elevación mediante elevatr. Luego los recortamos, les reducimos la resolución y el tamaño, y los transponemos (los rasters y los arreglos numéricos vienen con orientaciones distintas en R).
# convertir en array
test_c_arr <- as.array(grRastst)
# descargar elevación
Rch_dem <- elevatr::get_elev_raster(grRastst, z = 9)
Rch_dem <- raster::crop(Rch_dem, grRastst)
plot(Rch_dem)
# reducir tamaño
areademMatrix <- rayshader::resize_matrix(rayshader::raster_to_matrix(Rch_dem), scale = 0.7)
# transponer
demFinal <- (areademMatrix)Así se ve el DEM.
Ahora calculamos el sombreado y combinamos el arreglo con los datos de elevación.
# Sombreado
ambient_layer <- ambient_shade(demFinal, zscale = 10, multicore = TRUE, maxsearch = 200)
ray_layer <- ray_shade(demFinal, zscale = 20, multicore = TRUE)
# 3D
(test_c_arr / 255) %>%
add_shadow(ray_layer, 0.3) %>%
add_shadow(ambient_layer, 0) %>%
plot_3d(demFinal,
zscale = 150, theta = -50, phi = 28, zoom = 0.3,
windowsize = c(1200, 800), solid = FALSE,
background = "#3d314a"
)Así se ve la ventana RGL que genera plot3d. Ajusté algunos de los valores de la cámara. Sale una vista que conocemos bien los chilangos.
También se ve bien la vista superior que resulta de plot_map()
# en 2D
(test_c_arr / 255) %>%
add_shadow(ray_layer) %>%
add_shadow(ambient_layer, 0) %>%
plot_map()
Finalmente, podemos añadir elementos con rayrender antes de exportar imágenes de alta resolución con render_highquality(). Me gustó esta alternativa para darle más ambiente a la escena con el argumento light = FALSE y una esfera de luz de color colocada en donde nos guste más.
Turn off the lights and add your own🤓: #rstats #rayshader #rayrender
— Tyler Morgan-Wall (@tylermorganwall) October 27, 2019
render_highquality(light = FALSE, scene_elements = sphere(y = 100, radius = 10, material = diffuse(lightintensity = 250, implicit_sample = TRUE))) pic.twitter.com/g1QDEphZGh
# render
render_highquality(
light = FALSE,
scene_elements = rayrender::sphere(
z = 450, y = 100, x = 200, radius = 36,
material = rayrender::light(
intensity = 150,
spotlight_width = 20,
color = "#ff773d"
)
)
)Aquí hay algunos renders que saqué con diferentes opciones de cámara e iluminación.
Los resultados me gustaron bastante, y con todo ésto ya hasta se pueden hacer animaciones. Esta manera de georeferenciar imágenes también sirve para mapas que hacemos con R base, siempre y cuando el panel ocupe todo el espacio de la imagen. Espero que esta guía les sirva. Quedo atento a cualquier duda o comentario. Salu-2