Creando Mapas de República Dominicana en R

Posted by Johan Rosa on Wednesday, October 23, 2019

TOC

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)

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()

Resaltar áreas en base a algún criterio

plot(map_region)
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")