En este documento se presentan algunos ejemplos de cartografía digital con R. Con el objetivo de entender principalmente el manejo de Sistemas de Información Geográfica (SIG) y aprender a visualizar datos e información que requiera de una localización específica. Simplemente presentando algunos objetivos y los pasos técnicos necesarios para llevarlos a cabo.
La información de la estructura de los mapas fue principalmente tomado por el Instituto Mora.

Sistema de Información Geográfica

El INEGI tiene puestos a disposición del público shapefiles con distintos niveles de información. Sin considerar los históricos esta información se regular por el Marco Geoestadístico Censal, cuya última versión es de 2019. Incluye las siguientes divisiones y subdivisiones del territorio:

Table: Subdivisiones territoriales básicas del INEGI.
Código	División	Unidades territoriales
AGEE	Estatales	32
AGEM	Municipales	2 465
AGEB	Censal	62 732
PLUR	Localides	295 359

Marco Geoestadístico

El Marco Geoestadístico es un sistema único y de carácter nacional diseñado por el Instituto Nacional de Estadística y Geografía (INEGI), el cual presenta la división del territorio nacional en diferentes niveles de desagregación para referir geográficamente la información estadística de los censos y encuestas institucionales y de las Unidades del Estado, que se integra al Sistema Nacional de Información Estadística y Geográfica (SNIEG).

Conformación del Marco Geoestadístico Nacional (MGN)

El Marco Geoestadístico Nacional está conformado por áreas geoestadísticas divida en tres áreas de desagregación:

Áreas geoestadísticas estatales (AGEE)
Áreas geoestadísticas municipales (AGEM)
Áreas geoestadísticas básicas (AGEB)
\(~\) Área geoestadística básica rural
\(~\) Área geoestadística básica urbana

Aspectos fundamentales

Dadas las diferencias de densidad de población y uso del suelo se consideran necesario distinguir dos tipos de AGEB: urbanas y rurales.

Las AGEB urbanas delimitan una parte o el total de una localidad de 2 500 habitantes o más, o bien, una cabecera municipal, independientemente de su número de pobladores, en conjuntos que generalmente van de 5 a 50 manzanas.

Las AGEB rurales enmarcan una superficie de aproximadamente 10 000 hectáreas, cuyo uso del suelo es predominantemente agropecuario y en ellas se encuentran distribuidas las localidades menores a 2 500 habitantes, que, para fines operativos, se han denominado como localidades rurales.

Al interior de estos niveles de desagregación, se encuentran las áreas geográficas que contienen las unidades mínimas de observación del Censo de Población y Vivienda 2010 (habitantes y viviendas) las cuales son:

Localidad
Manzana

El archivo .shp INEGI Marco Geoestadístico,2018, que se pueden descargar directamente o se puede correr el algoritmo de descarga.

#Algoritmo de descarga de los datos .zip
url<-"http://internet.contenidos.inegi.org.mx/contenidos/productos/prod_serv/contenidos/espanol/bvinegi/productos/geografia/marcogeo/889463674658_s.zip"
if(!file.exists("areas_geoestadisticas_estatales.zip"){
  download.file(url,"areas_geoestadisticas_estatales.zip",mode="wb")
}
if(!file.exists("conjunto_de_datos/areas_geoestadisticas_basicas_rurales.dbf")){
  unzip("areas_geoestadisticas_estatales.zip")
}

Mapas de México con R

General

Se lee el archivo shp que contiene la división política de México.

La función readOGR del paquete rgdal, que extrae automáticamente esta información es utilizada por otros paquetes SIG de código abierto como QGIS y permite a R manejar una gama más amplia de formatos de datos espaciales.

shape_estados <- readOGR(dsn ="MGN/MGN Junio 2016/conjunto_de_datos", 
                          layer = "areas_geoestadisticas_estatales",
                          #  encoding = "UTF-8",
                              use_iconv = TRUE)

class(shape_estados)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

proj4stringrepresenta el sistema de referencia de coordenadas utilizado en los datos.

proj4string(shape_estados)

## [1] "+proj=lcc +lat_0=12 +lon_0=-102 +lat_1=17.5 +lat_2=29.5 +x_0=2500000 +y_0=0 +ellps=GRS80 +units=m +no_defs"

Mapa simple de México

ggplot(data=shape_estados)+
  geom_polygon(aes(x=long, 
                    y=lat, 
                      group=group)) + 
    labs(title = "Mapa simple de México")

División política de los estados

ggplot(data=shape_estados) +  
  geom_polygon(aes(x=long, y=lat, group=group), 
                  fill="white",
                    color="black") +
      theme_minimal() +
        theme(axis.text = element_blank()) +
   labs(title = "Mapa simple de México")

Mapas coropléticos de México con R.

Un mapa de coropletas muestra áreas o regiones geográficas divididas que están coloreadas en relación con una variable numérica o categórica.

Para este caso, se toma como ejemplo el índice de marginación 2015 elaborado por el Consejo Nacional de Población (CONAPO), donde este tipo de índice enfatiza la cuestión territorial, la población que vive territorios marginados.

Índice de Marginación 2015 por Estado

La base de datos del índice de marginación por estado se encuentra disponible en la página de datos abiertos.

IM_2015<-read.csv("http://www.conapo.gob.mx/work/models/CONAPO/Marginacion/Datos_Abiertos/Entidad_Federativa/Marginacion_por_entidad_2015.csv", header=TRUE, sep=",") %>%
  mutate (POBTOT2015=gsub(" ","", POBTOT2015)) %>%
        mutate(POBTOT2015=as.integer(POBTOT2015),
               CVE_ENT=as.character(.$CVE_ENT),
               CVE_ENT=ifelse(nchar(CVE_ENT)==1,
                                  paste0("0",CVE_ENT),
                                  CVE_ENT))%>%
            filter(NOM_ENT!="Nacional")%>%
              select(.,1:15) %>%
                mutate(id=as.character(seq(0,to=(nrow(.)-1))))

La función fortify() del paquete ggplot2, permite convertir un objeto (ej.spdf) a un data.frame. Al hacer este tipo de conversión permite unir los datos de Índice de marginación que es un data.frame al SpatialPolyDataFame.

Se crea una variable capas_estados, en la cual se convirtió el spdf a un data.frame y se agrego la variable grado de marginación 2015 GME, utilizando la función inner_join.

capas_estados<-inner_join(fortify(shape_estados,by="CVE_ENT"),IM_2015%>%select(id,GME),by="id")
names(capas_estados)

## [1] "long"  "lat"   "order" "hole"  "piece" "id"    "group" "GME"

geom_polygon() dibuja líneas entre puntos y “los cierra” (es decir, dibuja una línea desde el último punto hasta el primer punto).

values<-c('tomato1', 'sandybrown','wheat3','darkolivegreen3','green3')

capas_estados %>% 
   mutate(GME=fct_relevel(.$GME,"Muy alto","Alto","Medio","Bajo","Muy bajo")) %>%
     ggplot() +  
       geom_polygon(aes(x=long, 
                         y=lat, 
                           group=group, #El argumento group=group arma grupos tanto para para los polígonos. 
                             color="black",
                               fill=GME)) + 
          theme_void() + 
            theme(plot.title = element_text(size=22),
                  legend.key.size = unit(0.5, "cm"),
                  legend.position = c(0.8, 0.7)) + 
              scale_fill_manual(values=values) + 
                scale_color_manual(values=c("black")) +
                  guides(color = FALSE) +
        labs(title = "Índice de marginación a nivel estatal",
             fill = "GM",
             caption = "Fuente: Encuesta Intercensal 2015 (INEGI), \n Índice de marginación por estado 2015 (CONAPO).")

Múltiples estructuras

Usando `geom_polygon()`

El paquete de ggplot2 tiene un elemento geométrico específico para generar mapas. Usa el primitivo geom_polygon(), como lo hicimos en el ejercicio anterior, pero permite especificar directamente el mapa y automatiza el manejo de múltiples estructuras de datos.

capas_estados %>% 
   mutate(GME=fct_relevel(.$GME,"Muy alto","Alto","Medio","Bajo","Muy bajo")) %>%
     ggplot() +  
       geom_polygon(aes(x=long, 
                         y=lat, 
                           group=group, #El argumento group=group arma grupos tanto para para los polígonos. 
                             color="black",
                               fill=GME)) + 
          theme_void() + 
            theme(plot.title = element_text(size=22),
                  legend.key.size = unit(0.5, "cm"),
                  legend.position = c(0.8, 0.7)) + 
              scale_fill_viridis_d() + 
                scale_color_manual(values=c("#BDBDBD")) +
                  guides(color = FALSE) +
      labs(title = "Índice de marginación a nivel estatal",
           subtitle = "Utilizando geom_polygon",
           fill = "GM",
           caption = "Fuente: Encuesta Intercensal 2015 (INEGI), \n Índice de marginación por estado 2015 (CONAPO).")

Usando `geom_map()`

geom_map generalmente será mucho más rápido que geom_polygon porque puede hacer la fusión de una manera más eficiente, uniendo la base de datos con el shape mediante el argumento map_id.

IM_2015 %>%
  mutate(GME=fct_relevel(.$GME,"Muy alto","Alto","Medio","Bajo","Muy bajo")) %>%
    ggplot(aes(map_id = id)) + 
      geom_map(aes(fill = GME,
                   color="GME"),
                    map = capas_estados)  +
        geom_errorbarh(aes(x=1.5e6, 
                           xmin=1.5e6-1e5,
                           xmax=1.5e6+1e5, y=1e6), 
                       height=5e4) +
          expand_limits(x = capas_estados$long, 
                        y = capas_estados$lat) +
           theme_void() + 
              theme(plot.title = element_text(size=22),
                    legend.key.size = unit(0.5, "cm"),
                    legend.position = c(0.8, 0.7)) + 
                 scale_fill_viridis_d(option="B",begin=0.2,end=0.9,direction = -1) + 
                   scale_color_manual(values=c("#BDBDBD")) +
                     guides(color = FALSE) + 
         annotate("text", x= 1.5e6, y=1.0e6+8e4, label="200km", size=2) +
            labs(title = "Índice de marginación a nivel estatal",
                 subtitle = "Utilizando geom_map",
                 fill = "GM",
                 caption = "Fuente: Encuesta Intercensal 2015 (INEGI), \n Índice de marginación por estado 2015 (CONAPO).")

Usando `tm_shape()`

La ventaja de usar tmap es que comparte una lógica con ggplot2. Por lo tanto, primero debemos crear un objeto tmap (tm_shape()) seguido de una capa temática (tm_polygons()), y luego se deben especificar parámetros adicionales (por ejemplo, tm_layout()). Los elementos siguientes los agregamos usando un +.

Para el mapa necesitamos un shapefile con las delimitaciones de los estados y un archivo .csv con valores del grado de marginación 2015, ambos almacenados en el ambiente [shape_estados y IM_2015]. Utilizando la función merge(), ayuda a fusionar un SpatialPolygonsDataFrame y un data.frame, lo que genera un spdf.
Anteriormente se transformó el shape [shape_estados]spdf en un data.frame, utilizando la función fortify().

require(tmap)
require(tmaptools)

# Fusionamos un df y un spdf
capas_estados<-merge(shape_estados,
                     IM_2015%>% 
                          mutate(GME=fct_relevel(.$GME,"Muy alto","Alto","Medio","Bajo","Muy bajo")),
                      by = "CVE_ENT")

tm_shape(capas_estados) + 
  tm_polygons("GME",
                popup.vars=c("GME:" = "GME",
                             "Population:" = "POBTOT2015"),
                  id = "CVE_ENT",
                    alpha = 1, # transparency, number of classes and palette
                      n = 5,
                        title = "GM",
                          palette = hcl.colors(18, palette = "Inferno",)[3:9]) + 
  tm_borders("grey20") +
    tm_layout(main.title = "Ìndice de marginación a nivel estatal", 
              frame=FALSE) + 
      tm_credits("Fuente: Encuesta Intercensal 2015 (INEGI), \nÍndice de marginación por estado 2015 (CONAPO).",
                 position=c("left", "bottom"))

Estructura de mapas

El sistema de coordenadas que utiliza INEGI, también definido en el el Marco Geoestadístico Censal 2016. La definición técnica es MEXICO_ITRF_2008_LCC, con proyección Cónica Conforme de Lambert(CCL) y datum es ITRF2008. Se trata de un sistema de coordenadas rectangular, cuya unidad de medida es el metro. Esta cuadricula cubre el espacio entre los paralelos 17.5N y 29.5N con el meridiano 102 como centro. En términos prácticos esto tiene algunas implicaciones:

La cartografía que generemos a partir de shapefiles del INEGI se expresa directamente en una proyección ad hoc para México. Por lo tanto:
\(~\) - Podemos esperar mapas con pocas deformaciones.
\(~\) - Al utilizar unidades rectas para expresar la distancia (y no grados, minutos, segundos, unidades de arco) es fácil interpretar la escala de los gráficos. Esta es razonablemente homogénea.
\(~\) - La mayor parte de la información geográfica que producen las instituciones gubernamentales de México (INEGI, CONAPO, CONABIO, etc.) utiliza este mismo marco de referencia, por lo que podemos combinar estas capas sin mayores inconvenientes. El INE es un caso aparte, su sistema de información geográfica es diferente.
Sin embargo otras capas de información geográfica pudiéramos aprovechar y que se produce fuera de la órbita de estas instituciones suele estar expresada en otro sistema de coordenadas. Por ejemplo, en unidades hexadecimales de arco, es decir, latitud y longitud medidas en grados, minutos y segundos. Por lo tanto:
\(~\) - No podemos unir esas capas sin una transformación previa, de lo contrario la ubicación de los polígonos no va a ser compatible.
\(~\) - Debemos llevar la información no INEGI al sistema de INEGI o viceversa.
\(~\) - INEGI dispone de una app en línea para convertir información en lat-long \(\rightarrow\) LCC en http://www.inegi.org.mx/geo/contenidos/geodesia/traninv.aspx. Hasta ahora no hay una opción para convertir una base de puntos datos completa y mucho menos polígonos. Debemos transformar punto por punto.
\(~\) - No encontré una función o librería en R que nos permita hacer este traspaso de manera simple, aunque no he sido exhaustivo en la búsqueda.
\(~\) - http://www.gadm.org/ dispone de shapefiles y objetos geoespaciales SPDF de R (.rds) de México, subdivididos hasta nivel municipal con un sistema de coordenadas geográficas. Se trata de archivos pequeños que se cargan directamente en R sin necesidad de importación. Sin embargo no dispone de unidades territoriales menores al municipio, red caminera, áreas y recursos naturales, etc. y no es fuente oficial. Tampoco podemos estar completamente seguros de su actualización. Úselo bajo su propia responsabilidad.
\(~\) - La cartografía de GoogleMaps está disponible a través de una API y R tiene varios clientes para esa API. El problema señalado de los diferentes sistemas de coordenadas dificulta su aprovechamiento.
- R tiene la librería mapproj:: para proyección de mapas y ggplot nos permite utilizarla fácilmente con coord_map(). Sin embargo, esta librería sólo admite coordenadas geográficas, por lo que no podemos utilizarla con datos crudos del INEGI.

Mapas coropléticos de México con R.

Los mapas coropléticos representan polígonos de unidades territoriales a los que se asigna un color de acuerdo al valor de una variable pertinente. Este color puede representar tanto una escala continua-en la que se genera un degradado continuo de color- o una escala discreta -en la que cada categoría de esa variable tiene un color asignado.

El procedimiento en R usando ggplot2:: es el siguiente:

1. Importar los polígonos con los contornos de las unidades territoriales. Esto producirá un objeto de la clase SPDF o Spatial Polygon Data Frame.

\(~\) - Los objetos SPDF NO son directamente un data.frame, en sentido estricto pertenecen a la clase list. Sin embargo se transforman fácilmente a data.frame.

\(~\) - Fueron incorporados en R a partir de la librería sp::, un conjunto de clases y métodos para datos espaciales.
\(~\:\) -sp:: es la librería primitiva para información geográfica en R, prácticamente todas las demás depende de un modo u otro de las especificaciones, funciones y métodos de sp::.

2. Unir los datos estadísticos que controlarán el color con la estructura de datos poligonal.
\(~\) - A través de un join_inner o con el método simplificado de geom_map().

3.Generar el mapa con una llamada a ggplot() y los elementos geométricos geom_polygon() o geom_map() del paquete ggplot2::.

Extensión geográfica

Se utiliza la función bbox del paquete sp:: donde este trata de identificar la extensión geográfica de los mosaicos de imágenes que necesitamos del shape.

b_box<-bbox(shape_estados)

	min	max
x	911292.0	4082997
y	319149.1	2349615

Librerías

Librerías que se usaron en el trabajo

##  [1] "tmaptools"  "tmap"       "forcats"    "ggplot2"    "rgdal"     
##  [6] "sp"         "foreign"    "dplyr"      "knitr"      "kableExtra"

Referencias

Interactive choropleth maps with R and tmap (part I) | Marcin Stępniak. (2019). Retrieved November 2, 2020, from https://marcinstepniak.eu/post/interactive-choropleth-maps-with-r-and-tmap-part-i/

Marco Geoestadístico. (2020). Retrieved November 1, 2020, from https://www.inegi.org.mx/temas/mg/

Paladino, M. (2017). Mapas con ggplot. Retrieved November 1, 2020, from https://www.institutomora.edu.mx/testU/SitePages/martinpaladino/Mapas_con_R.html

tmap: get started! (2020). Retrieved November 3, 2020, from https://cran.r-project.org/web/packages/tmap/vignettes/tmap-getstarted.html

This work by Diana Villasana Ocampo is licensed under a Creative Commons Attribution 4.0 International License.

Mapa por estados de México

Diana Villasana Ocampo

Sistema de Información Geográfica

Marco Geoestadístico

Mapas de México con R

General

Mapa simple de México

División política de los estados

Mapas coropléticos de México con R.

Índice de Marginación 2015 por Estado

Múltiples estructuras

Usando geom_polygon()

Usando geom_map()

Usando tm_shape()

Estructura de mapas

Mapas coropléticos de México con R.

Extensión geográfica

Librerías

Referencias

Usando `geom_polygon()`

Usando `geom_map()`

Usando `tm_shape()`