Una manera muy interesante de visualizar los datos es mediante mapas, vestir la cartografía con información. Los mapas constituyen una herramienta muy poderosa para comunicar, ya que llaman la atención de los lectores y suelen ser interpretados casi intuitivamente.
Los software tradicionales para hacer este tipo de cosas son Qgis, ArcGIS y otros, pero si estás aprendiendo R resultaría muy útil dominar algunos paquetes que no tienen nada que envidiarles a las herramientas mencionadas.
Paqutes
library(tmap) # graficar objetos de poligonos espaciales
library(rgdal) # lectura de objetos "spatialpolygonsdataframe"
library(tidyverse) # visualización y manipulación de datos
#library(xlsx) # leer archivos de excel
library(kableExtra) # para dar formato a las tablas
library(leaflet) #mapas interactivos
library(raster) # descargar shape files (también otras cosas)
library(gridExtra) # poner variaor gráficos juntos
library(janitor) # limpiar datos
library(sf) # Para manejar un tipo de data espacial
library(viridis) # Mapas
objetos tipo SpatialPolygonsDataFrame
En R, los objetos SpatialPolygonsDataFrame
son una especie de lista que contiene los objetos necesarios para hacer mapas. Tradicionalmente a estos objetos se les denomina shapefiles y tienen un componente con los polígonos de las entidades geográficas, un componente que almacena la información referida al sistema de coordenadas geográficas de los polígonos, un data frame con una fila para cada entidad geográfica del shapefile, y un metadato sobre la proyección geográfica con la que se construyó el archivo.
Los shapefiles de la cartografía dominicana se pueden encontrar en la página web de la Oficina Nacional de Estadística. Hay archivos con los polígonos a diferentes niveles de desagregación, desde regiones y provincias, hasta llegar a barrios, el nivel mínimo. Con el siguiente enlace pueden acceder a ellos. Repositorio personal o web ONE.
Importar, visualizar y manipular shapefiles
Ya que tenemos los archivos .shp
de la cartografía nacional, podemos importar los correspondientes a regiones de planificación. Para esto usaremos la función readOGR()
.
map_region <- readOGR(dsn = "/cloud/project/mapas_rd/shape_files/REGCenso2010.shp")
map_provincia <- readOGR(dsn = "/cloud/project/mapas_rd/shape_files/PROVCenso2010.shp")
Una vez importado el objeto y nombrado en el ambiente de trabajo como map_provincia
, mostrar el mapa resulta muy sencillo, ya que la función básica plot
de R
puede manejar este tipo de archivos.
plot(map_region)
## Warning in wkt(obj): CRS object has no comment
Accediendo al componente data
El componente data de los SpatialPolygonsDataFrame
es el elemento que con más frecuencia tendremos que manipular. Esto porque por lo general los shapefiles no traen información sobre las áreas geográficas más que aquellas sirven para identificar las mismas.
map_region
simplemente tiene el código de la región y el nombre, pero más adelante le agregaremos otros datos.
map_region@data %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condense"), full_width = T)
REG | region | poblacion | pextrema | pgeneral |
---|---|---|---|---|
01 | Cibao Norte | 1600820 | 2.392323 | 18.79252 |
02 | Cibao Sur | 733379 | 2.100254 | 14.22712 |
03 | Cibao Nordeste | 641259 | 0.592932 | 16.84867 |
04 | Cibao Noroeste | 413673 | 2.985500 | 30.54892 |
05 | Valdesia | 1096613 | 4.356706 | 26.42024 |
06 | Enriquillo | 381158 | 6.860147 | 45.34002 |
07 | El Valle | 287657 | 7.803186 | 42.42642 |
08 | Yuma | 698277 | 2.614538 | 19.90459 |
09 | Higuamo | 578478 | 3.983395 | 25.31570 |
10 | Ozama | 3834835 | 2.323483 | 21.55243 |
Para agregar información al componente data del shapefile, vamos a descargar unos cuadros que publica la ONE con las proyecciones de población y las tasas de pobreza por regiones. Estos archivos fueron creados con fines de impresión y aunque están en excel contienen headers y subheaders en celdas unificadas, por lo que se necesitará limpiar la data. Los pasos para hacer esto los compartiré, con fines de garantizar la reproducibilidad de todo el contenido, pero no daré muchos detalles porque el documento es sobre visualización cartográfica más que manipulación de datos con R.
# Descargar archivos ---------------------------------------------------------------
url_poblacion <- "https://www.one.gob.do/Multimedia/Download?ObjId=29649"
url_pobreza <- "https://www.one.gob.do/Multimedia/Download?ObjId=88791"
download.file(url_poblacion, "mapas_rd/proy-poblacion.xlsx")
download.file(url_pobreza, "mapas_rd/tasa_pobreza.xlsx")
# Importar archivos ---------------------------------------------------------------
poblacion <- read.xlsx(
"/cloud/project/mapas_rd/proy-poblacion.xlsx",
sheetIndex = 1,
startRow = 9,
header = F
)
pobreza <- read.xlsx(
"mapas_rd/tasa_pobreza.xlsx",
sheetIndex = 1,
startRow = 7,
endRow = 26,
header = F
)
# Adecuar para hacer join con la data del shapefile --------------------------------
# crear header en base a los archivos originales, pero adecuados para ser usados
# en r. En otro post del blog les mostraré como adecuar estos header de otra manera
header_poblacion <- c(
"region",
paste(
rep(c("total", "hombre", "mujer"), times = 31),
rep(2000:2030, each = 3), sep = "_")
)
header_pobreza <- c(
"year",
paste(
rep(
c("Cibao Norte", "Cibao Sur", "Cibao Noreste", "Cibao Noroeste", "Valdesia",
"Enriquillo", "El Valle", "Yuma", "Higuamo", "Ozama", "Pais"),
times = 2),
rep(c("pgeneral", "pextrema"), each = 11), sep = "_")
)
# Agregando los nuevos headers
names(pobreza) <- header_pobreza
names(poblacion) <- header_poblacion
Vamos a utilizar los datos de población y tasas de pobreza para el 2018 para agregar datos al shapefile
Proyecciones de población - 2018
poblacion2018 <-
poblacion %>%
gather(sexo, poblacion, -region) %>%
separate(col = sexo, into = c("sexo", "year")) %>%
filter(!is.na(poblacion)) %>%
filter(!region == "Total país") %>%
filter(str_detect(region, "Región"),
sexo == "total",
year == "2018") %>%
dplyr::select(-sexo, -year) %>%
mutate(region = str_remove(region, "^Región ") %>%
str_trim(),
region = recode(
region,
"Metropolitana" = "Ozama",
"Del Valle" = "El Valle"))
# Mostrar la tabla
poblacion2018 %>%
mutate(
poblacion = format(poblacion, big.mark = ",")) %>%
arrange(region) %>%
kable() %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condense"),
full_width = T)
region | poblacion |
---|---|
Cibao Nordeste | 641,259 |
Cibao Noroeste | 413,673 |
Cibao Norte | 1,600,820 |
Cibao Sur | 733,379 |
El Valle | 287,657 |
Enriquillo | 381,158 |
Higuamo | 578,478 |
Ozama | 3,834,835 |
Valdesia | 1,096,613 |
Yuma | 698,277 |
pobreza general y extrema - 2018
pobreza2018 <-
pobreza %>%
filter(
!is.na(`Cibao Norte_pgeneral`),
year == "2018 ENCFT") %>%
gather(region, tasa_pobreza, -year) %>%
separate(
region, c("region", "tipo_pobreza"), "_") %>%
spread(tipo_pobreza, tasa_pobreza) %>%
filter(!region == "Pais") %>%
mutate(
region = recode(region, "Cibao Noreste" = "Cibao Nordeste")) %>%
dplyr::select(-year)
# Mostrar la tabla
pobreza2018 %>%
arrange(region) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condense"), full_width = TRUE)
region | pextrema | pgeneral |
---|---|---|
Cibao Nordeste | 0.592932 | 16.84867 |
Cibao Noroeste | 2.985500 | 30.54892 |
Cibao Norte | 2.392323 | 18.79252 |
Cibao Sur | 2.100254 | 14.22712 |
El Valle | 7.803186 | 42.42642 |
Enriquillo | 6.860147 | 45.34002 |
Higuamo | 3.983395 | 25.31570 |
Ozama | 2.323483 | 21.55243 |
Valdesia | 4.356706 | 26.42024 |
Yuma | 2.614538 | 19.90459 |
Para agregar información a la data de los shapefiles simplemente haremos un join por región.
# Join de los 3 dataframes
map_region@data <- map_region@data %>%
rename(region = TOPONIMIA) %>%
mutate(
region = str_to_title(region) %>%
str_remove("^Región ") %>%
#str_trim() %>%
recode("Ozama O Metropolitana" = "Ozama")) %>%
left_join(poblacion2018) %>%
left_join(pobreza2018)
# Mostrar la tabla
map_region@data %>%
mutate_at(.vars = c("poblacion", "pextrema", "pgeneral"), round, digits = 2)%>%
mutate_at(.vars = c("poblacion", "pextrema", "pgeneral"), format, big.mark = ",") %>%
arrange(region) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "hover", "condense"), full_width = TRUE)
REG | region | poblacion | pextrema | pgeneral |
---|---|---|---|---|
03 | Cibao Nordeste | 641,259 | 0.59 | 16.85 |
04 | Cibao Noroeste | 413,673 | 2.99 | 30.55 |
01 | Cibao Norte | 1,600,820 | 2.39 | 18.79 |
02 | Cibao Sur | 733,379 | 2.10 | 14.23 |
07 | El Valle | 287,657 | 7.80 | 42.43 |
06 | Enriquillo | 381,158 | 6.86 | 45.34 |
09 | Higuamo | 578,478 | 3.98 | 25.32 |
10 | Ozama | 3,834,835 | 2.32 | 21.55 |
05 | Valdesia | 1,096,613 | 4.36 | 26.42 |
08 | Yuma | 698,277 | 2.61 | 19.90 |
Mapas estáticos
Primero vamos a explorar un poco más los gráficos con la función base plot()
. Para esto vamos a hacer dos cosas, filtrar áreas geográficas y graficarlas (Aunque tenemos el mapa nacional con todas la regiones, no siempre querríamos graficar todas) y Colorear en el mapa regiones específicas.
La lógica de filtrar los shape files es muy similar a la de filtrar un dataframe, aquí mostraremos como seleccionar y graficar las áreas de la región norte o Cibao.
Para resaltar áreas geográficas específicas simplemente se utilizar una condición lógica similar a la del mapa anterior, pero antes se grafica el mapa general y se le agrega un layer con los colores de las zonas a resaltar. A continuación muestro el código para realizar ambas acciones.
Filtrar el shapefiles
map_region[str_detect(map_region$region, "Cibao"),] %>%
plot()
## Warning in wkt(obj): CRS object has no comment
Resaltar áreas en base a algún criterio
plot(map_region)
## Warning in wkt(obj): CRS object has no comment
plot(map_region[map_region$pgeneral > mean(map_region$pgeneral),], col = "turquoise", add = T)
Otro tipo de mapa estático que se puede realizar es colocando colores a los polígonos dependiendo de intervalos de alguna variable, de manera que si se usa la escala de color adecuada permite identificar las distribución de esta en las diferentes zonas. Para hacer este tipo de mapas recurriremos al paquete tmap
, el cuál es un paquete basado el la gramática de gráficos en capas en la que se basa ggplot2
.
Mapa plano con tmap
tm_shape(map_region) +
tm_polygons()
Colores en base a una variable
tm_shape(map_region) +
tm_polygons("poblacion")
Ajustar los intervalos
tm_shape(map_region) +
tm_polygons("poblacion",
palette = "Blues", #"reds", "Puerples",
breaks = c(0, 500000, 1000000, 2000000,4000000))
Mapas interactivos
Los mapas interactivos se crean por lo general en formato HTML y permiten que el usuario cambie aspectos de los mismos, o bien, ellos cambian automáticamente dependiendo de algunas acciones. Entre las cosas más comunes que admiten este tipo de objetos están la posibilidad de modificar la escala del gráfico (Zoom), desplegar leyendas y etiquetas al desplazar el puntero por el gráfico y la aparición o desaparición de elementos en la medida que se modifica la escala del gráfico.
La librería que utilizo por excelencia para hacer este tipo de mapas es leaflet
la cual admite shape files o crea gráficos en base a openstreetmap.
Hagamos dos mapas con leaflet
, primero uno solo con los layers de openstreetmap y luego otro con los polígonos que venimos utilizando desde el principio de la publicación.
Lo primero que hay que hacer para utilizar los shapefiles que descargamos de la ONE es cambiar el sistema de coordenadas geográficas con el que viene por defecto. La mayoría de paquetes para hacer mapas usan principalmente la coordenadas basadas en grados, minutos y segundos; pero el objeto map_region
tiene una proyección distinta la UTM
A continuación coloco el comando para transformar la proyección del mapa. Como ven, este es otro de los elementos que le podemos modificar a este tipo de objetos.
Para aprender un poco más sobre los sistemas de referencia y coordenadas puede consultar este enlace
map_region <- spTransform(map_region, CRS("+proj=longlat +datum=WGS84"))
map_provincia <- spTransform(map_provincia, CRS("+proj=longlat +datum=WGS84"))
Mapa de las regiones
leaflet() %>%
addPolygons(
data = map_region,
color = "gray")