10  Librería leaflet

10.1 Introducción

Leaflet es una de las bibliotecas de JavaScript de código abierto más populares para crear mapas interactivos. Instituciones y empresas de índole muy diversa, como la ONG Acción contra el hambre o Flickr, así como especialistas en SIG como OpenStreetMap, utilizan Leaflet para mostrar información geográfica en sus respectivos sitios web.

10.2 Características

  • Panorámica/zoom interactivo.
  • Componer mapas utilizando combinaciones arbitrarias de:
    • Mosaicos de mapa
    • Marcadores
    • Polígonos
    • Líneas
    • Ventanas emergentes
    • GeoJSON
  • Crear mapas directamente desde la consola R o RStudio.
  • Incrustar mapas en documentos knitr/R Markdown y aplicaciones Shiny.
  • Representar fácilmente objetos espaciales de los paquetes sp o sf, o marcos de datos con columnas de latitud/longitud.
  • Utilizar los límites del mapa y los eventos del mouse para impulsar la lógica Shiny.
  • Mostrar mapas en proyecciones Mercator no esféricas.

10.3 Instalación

Para instalar el paquete leaflet(Cheng et al., 2024), usamos la función habitual install.packages():

install.packages("leaflet")

Si lo deseamos, podemos instalar la versión en desarrollo (no es necesario para esta práctica):

devtools::install_github("rstudio/leaflet")

Una vez instalada la librería, ya se puede cargar y usar desde R, en documentos R Markdown o Quarto (como este) o en aplicaciones Shiny como esta.

library(leaflet)

10.4 Funciones básicas

Antes de comenzar, cargamos algunas de las librerías habituales:

library(tidyverse)
library(sf)
library(classInt)

10.4.1 Mapa centrado en un lugar

Puesto que lo que queremos es representar algo en un mapa, lo mínimo que necesitamos es cargar el mapa y “pintar” algo en él.

Para cargar el mapa usamos la función addTiles().

Podemos pintar un punto, por ejemplo la ubicación de este edificio. Pero…

Pregunta

¿Cómo conseguimos las coordenadas de los marcadores que queremos mostrar?

leaflet() %>%
  addTiles() %>%
  addMarkers(lat=39.47776,
             lng=-0.34374)

10.4.2 Markers

Si queremos que al hacer clic sobre el icono aparezca el texto deseado debemos incluir en la función addMarkers() el parámetro popup:

leaflet() %>%
  addTiles() %>%
  addMarkers(lat=39.47776,
             lng=-0.34374,
             popup="Facultat d'Economia - UV")

Podemos usar circulos en lugar de iconos, usando la función addCircleMarkers():

leaflet() %>%
  addTiles() %>%
  addCircleMarkers(lat=39.47776,
                   lng=-0.34374,
                   popup="Facultat d'Economia - UV",
                   color="red",
                   radius = 8,
                   stroke = 2,
                   opacity = 0.5)

10.4.3 Zoom

Podemos centrar el mapa en otro lugar, usando el parámetro zoom de la función setView().

También debemos indicar las coordenadas que usaremos como centro del mapa.

El zoom por defecto es 15.

leaflet() %>%
  addTiles() %>%
  addMarkers(lat=39.47776,
             lng=-0.34374,
             popup="Facultat d'Economia - UV") %>%
  setView(lat=40.41428,
          lng=-3.72864,
          zoom=5)

10.4.4 Mapa base

Si queremos usar otro mapa distinto al cargado por defecto, usaremos la función addProviderTiles().

leaflet() %>%
  addProviderTiles("CartoDB.Positron") %>%
  addMarkers(lat=39.47776,
             lng=-0.34374,
             popup="Facultat d'Economia - UV")

Podemos seleccionar entre una amplia variedad de capas (Tiles). Hay más de 200 opciones disponibles. Aquí algunos ejemplos:

head(as.character(providers),20)
 [1] "OpenStreetMap"            "OpenStreetMap.Mapnik"    
 [3] "OpenStreetMap.DE"         "OpenStreetMap.CH"        
 [5] "OpenStreetMap.France"     "OpenStreetMap.HOT"       
 [7] "OpenStreetMap.BZH"        "MapTilesAPI"             
 [9] "MapTilesAPI.OSMEnglish"   "MapTilesAPI.OSMFrancais" 
[11] "MapTilesAPI.OSMEspagnol"  "OpenSeaMap"              
[13] "OPNVKarte"                "OpenTopoMap"             
[15] "OpenRailwayMap"           "OpenFireMap"             
[17] "SafeCast"                 "Stadia"                  
[19] "Stadia.AlidadeSmooth"     "Stadia.AlidadeSmoothDark"

En este enlace podemos visualizar todo el repertorio.

10.5 Funciones avanzadas

10.5.1 Puntos

Vamos a representar la ubicación de los centros educativos ubicados en la Comunitat Valenciana, utilizando el conjunto de datos que trabajamos (y mejoramos) en sesiones anteriores. El fichero con los datos, en formato GeoPackage, está disponible aquí.

data_edu <- sf::st_read("data/centres_educatius_bueno.gpkg", quiet=T)

Echemos un vistazo al dataset:

glimpse(data_edu)
Rows: 3,708
Columns: 34
$ codcen            <chr> "03000047", "03000060", "03000096", "030…
$ dlibre            <chr> "CEIP LA RAMBLA", "CENTRE PRIVAT LA MILA…
$ dgenerica_val     <chr> "COL·LEGI D'EDUCACIÓ INFANTIL I PRIMÀRIA…
$ despecifica       <chr> "LA RAMBLA", "LA MILAGROSA", "MIRADOR D'…
$ regimen           <chr> "PÚBLICO", "CONCERTADO", "PÚBLICO", "PÚB…
$ direccion         <chr> "AV. DE ALCOY, S/N", "CL. JUAN XXIII, 2"…
$ codpos            <chr> "03698", "03698", "03569", "03340", "033…
$ municipio_oficial <chr> "Agost", "Agost", "Aigües", "Albatera", …
$ provincia_val     <chr> "Alacant", "Alacant", "Alacant", "Alacan…
$ telef             <chr> "966908135", "965691143", "966908140", "…
$ fax               <chr> "0", "0", "0", "0", "0", "0", "0", "0", …
$ mail              <chr> "03000047@edu.gva.es", "cdlamilagrosa@co…
$ codauto           <chr> "10", "10", "10", "10", "10", "10", "10"…
$ iso2.ccaa.code    <chr> "ES-VC", "ES-VC", "ES-VC", "ES-VC", "ES-…
$ nuts1.code        <chr> "ES5", "ES5", "ES5", "ES5", "ES5", "ES5"…
$ nuts2.code        <chr> "ES52", "ES52", "ES52", "ES52", "ES52", …
$ ine.ccaa.name     <chr> "Comunitat Valenciana", "Comunitat Valen…
$ iso2.ccaa.name.es <chr> "Valenciana, Comunidad", "Valenciana, Co…
$ iso2.ccaa.name.ca <chr> "Valenciana, Comunitat", "Valenciana, Co…
$ iso2.ccaa.name.gl <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ iso2.ccaa.name.eu <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
$ nuts2.name        <chr> "Comunidad Valenciana", "Comunidad Valen…
$ cldr.ccaa.name.en <chr> "Valencian Community", "Valencian Commun…
$ cldr.ccaa.name.es <chr> "Comunidad Valenciana", "Comunidad Valen…
$ cldr.ccaa.name.ca <chr> "País Valencià", "País Valencià", "País …
$ cldr.ccaa.name.ga <chr> "Comunidade Valenciana", "Comunidade Val…
$ cldr.ccaa.name.eu <chr> "Valentziako Erkidegoa", "Valentziako Er…
$ ccaa.shortname.en <chr> "Valencian Community", "Valencian Commun…
$ ccaa.shortname.es <chr> "Comunidad Valenciana", "Comunidad Valen…
$ ccaa.shortname.ca <chr> "País Valencià", "País Valencià", "País …
$ ccaa.shortname.ga <chr> "Comunidade Valenciana", "Comunidade Val…
$ ccaa.shortname.eu <chr> "Valentziako Erkidegoa", "Valentziako Er…
$ nuts1.name        <chr> "ESTE", "ESTE", "ESTE", "ESTE", "ESTE", …
$ geom              <POINT [°]> POINT (-0.635608 38.43839), POINT …
unique(data_edu$regimen)
[1] "PÚBLICO"    "CONCERTADO" "PRIVADO"

Para poder representar datos espaciales (objetos sf) con leaflet, tenemos que modificar el CRS a WGS84:

data_edu <- sf::st_transform(data_edu, 4326)

Representamos los datos en el mapa:

data_edu %>%
  leaflet() %>% 
  addTiles() %>%
  addCircles()

Creamos la paleta de colores para diferenciar los centros según regimen:

tipos_regimen <- data_edu$regimen %>% unique()
colores <- c('#66c2a5','#fc8d62','#8da0cb') # https://colorbrewer2.org/
pal <- leaflet::colorFactor(colores, domain = tipos_regimen)

Probemos nuevamente:

data_edu %>%
  leaflet() %>% 
  addTiles() %>%
  addCircles(color = ~pal(regimen))

Identifiquemos cada punto con el nombre de cada centro:

m <- data_edu %>%
  leaflet() %>% 
  addTiles() %>%
  addCircles(color = ~pal(regimen),
             label = ~despecifica,
             popup = ~dgenerica_val,
             group = "Régimen")
m

Añadimos la leyenda:

m <- m %>%
  addLegend(data = data_edu,
            values = ~regimen,
            position = "bottomright",
            pal = pal,
            title = "Régimen",
            group = "Leyenda")
m

Añadimos control de capas:

m %>%
  addLayersControl(overlayGroups = c("Régimen","Leyenda"),
                   options = layersControlOptions(collapsed = T))

Sería interesante poder seleccionar qué centros visualizar (en función de la variable regimen)…

lista_regimen <- list()
for (i in 1:length(tipos_regimen)) {
  lista_regimen[[i]] <- data_edu %>% filter(regimen == tipos_regimen[i])
}
names(lista_regimen) <- tipos_regimen

map <- leaflet() %>% addTiles()
for (i in 1:length(lista_regimen)) {
  map <- map %>%
    addCircles(data = lista_regimen[[i]],
               color = ~pal(regimen),
               label = ~despecifica,
               popup = ~dgenerica_val,
               group = tipos_regimen[i])
}

map <- map %>%
  addLegend(data = data_edu,
            values = ~regimen,
            position = "bottomright",
            pal = pal,
            title = "Régimen",
            group = "Leyenda")

grupos <- c("Leyenda", tipos_regimen)

map2 <- map %>%
  addLayersControl(overlayGroups = grupos,
                   options = layersControlOptions(collapsed = T)) %>%
  hideGroup("Leyenda")
map2

Otra forma de hacerlo:

map3 <- map %>%
  addLayersControl(overlayGroups = "Leyenda",
                   baseGroups = tipos_regimen,
                   options = layersControlOptions(collapsed = T)) %>%
  hideGroup("Leyenda")
map3

También podemos visualizar los puntos de forma agrupada:

data_edu %>%
  leaflet() %>% 
  addProviderTiles("Esri.WorldImagery") %>%
  addCircleMarkers(popup = ~paste(despecifica, "<br>", regimen),
                   color = "red",
                   stroke = 2,
                   radius = 8,
                   clusterOptions = markerClusterOptions())

10.5.2 Polígonos

Vamos a representar la distribución de la renta de los hogares. Para ello, hemos escogido la Renta bruta media por hogar, para cada sección censal de la provincia de Valencia. Los microdatos (para el año 2021) se han obtenido del INE. También puedes descargar el fichero directamente de aquí

renta <- readxl::read_xlsx("data/renta_bruta_media_hogar_sscc_val_prov.xlsx")
glimpse(renta)
Rows: 1,823
Columns: 2
$ CUSEC <chr> "4600101001", "4600201001", "4600301001", "460040100…
$ renta <chr> "29986", "31872", "33446", "33916", "35777", "43522"…
renta <- renta %>% mutate(renta = as.numeric(renta))
Warning: There was 1 warning in `mutate()`.
ℹ In argument: `renta = as.numeric(renta)`.
Caused by warning:
! NAs introducidos por coerción
renta <- renta %>% filter(!is.na(renta))

La cartografía por secciones censales también la podemos obtener del INE (ver aquí).

Ejercicio
  1. Descarga y lee la cartografía del INE para el año 2021.
  2. Selecciona las secciones censales que pertenecen a la provincia de Valencia.
  3. Cambia el CRS a EPSG:4326.
  4. Guarda el dataset resultante como sscc_2021_CPRO_46.gpkg.

Unimos cartografía y datos de renta:

sscc_2021_val <- sf::st_read("data/sscc_2021_CPRO_46.gpkg", quiet = T)
renta_sf <- left_join(sscc_2021_val,renta, by="CUSEC")
Pregunta

¿Qué tendríamos que hacer si la variable de unión tuviera distinto nombre en cada dataset?

Creamos intervalos:

breaks <- classInt::classIntervals(var = renta_sf$renta, n=6, style = "quantile")
Warning in classInt::classIntervals(var = renta_sf$renta, n = 6,
style = "quantile"): var has missing values, omitted in finding
classes
breaks <- breaks$brks
breaks <- round(breaks, -2)

class <- base::findInterval(x=pull(renta_sf, renta), vec = breaks, all.inside = T)
renta_sf <- renta_sf %>% mutate(class = class)
labels <- tibble(class = 1:6,
               min = scales::number(breaks[1:6], big.mark = ".", decimal.mark = ","),
               max = scales::number(breaks[2:7], big.mark = ".", decimal.mark = ","),
               interval = paste0(min, " - ", max))

renta_sf <- inner_join(renta_sf, labels, by = "class")
renta_sf <- mutate(renta_sf, interval = factor(interval, levels = labels$interval))

En ocasiones, puede ser útil reducir la calidad de la cartografía (número de puntos que conforman los polígonos):

#install.packages("rmapshaper")
#library(rmapshaper)
renta_sf <- rmapshaper::ms_simplify(renta_sf)

Creamos la paleta de colores:

#install.packages("RColorBrewer")
#library(RColorBrewer)
#RColorBrewer::display.brewer.all()

palette <- RColorBrewer::brewer.pal(n=6, name = "YlGn")
colors <- tibble(clases = labels$interval, palette)
colors <- mutate(colors, clases = factor(clases, levels = labels$interval))
renta_sf <- mutate(renta_sf, clases = interval)
renta_sf <- inner_join(renta_sf, colors, by = "clases")

Representamos los datos en el mapa, con algunas funcionalidades extra que ofrece el paquete leaflet.extras(Gatscha et al., 2024):

#install.packages("leaflet.extras")
#library(leaflet.extras)

leaflet::leaflet() %>%
  leaflet::addTiles() %>%
  leaflet::addProviderTiles(provider = providers$CartoDB.Positron, group = "CartoDB") %>%
  leaflet::addPolygons(data = renta_sf,
              #stroke = T,
              color = "transparent",
              fillColor = renta_sf$palette,
              fillOpacity = 0.5,
              smoothFactor = 0.5,
              popup = paste0("Renta bruta media por hogar: ", renta_sf$renta, " €", "<br>", 
                             "Municipio: ", renta_sf$NMUN)) %>%
  leaflet::addLegend("topright",
            colors = colors$palette,
            labels = colors$clases,
            title = "Renta bruta media por hogar") %>%
  leaflet.extras::addResetMapButton() %>%
  leaflet.extras::addSearchOSM() %>%
  leaflet.extras::addControlGPS() %>%
  leaflet::addMiniMap(tiles = "Esri.WorldStreetMap", position = "bottomleft")

10.6 Recursos adicionales