Introducción

Con la reciente publicación de los resultados del censo 2020 de INEGI varios usuarios han estado presentando algunas gráficas en redes sociales con los principales resultados. Una de las visualizaciones más populares son los mapas coropléticos, especialmente a nivel municipal. Una ventaja de estos mapas es que nos permiten tener el panorama general de todo el país e incluso detectar a grandes rasgos algunos patrones regionales. Por ejemplo, acá @claudiodanielpc nos muestra el hacinamiento en los hogares:

Sin embargo, un inconveniente es que hay casos de regiones donde hay un gran número de municipios en una superficie pequeña lo que dificulta la lectura de la información. Un caso común en México es el área compuesta por los estados de Oaxaca, Puebla, Chiapas, Veracruz y Estado de México que contienen 1248 municipios, poco más de la mitad del total de 2469 municipios. Muchos mapas se vuelven difícil de ver en esta zona. Por ejemplo:

Para tratar de hacer más claras este tipo de visualizaciones podemos agregar uno o varios recuadros con acercamiento (zoom) a los resultados de estas zonas, tal como lo hace el mapa del ejemplo. Si bien eso no garantiza que toda la información se entenderá mejor, al menos puede ayudar en ciertas áreas. En esta entrada veremos cómo hacer varios tipos de recuadros para mapas en Python usando Geopandas y Matplotlib.

Empecemos importando las librerías necesarias.

import os
import sys

import geopandas as gpd
import matplotlib.pyplot as plt
import pandas as pd
import requests
import matplotlib_scalebar
from matplotlib_scalebar.scalebar import ScaleBar

print('Python', sys.version)
print(pd.__name__, pd.__version__)
print(gpd.__name__, gpd.__version__)
print(requests.__name__, requests.__version__)
print(plt.matplotlib.__name__, plt.matplotlib.__version__)
print(matplotlib_scalebar.__name__, matplotlib_scalebar.__version__)
Python 3.8.5 (default, Sep  3 2020, 21:29:08) [MSC v.1916 64 bit (AMD64)]
pandas 1.1.3
geopandas 0.8.1
requests 2.24.0
matplotlib 3.3.2
matplotlib_scalebar 0.7.2

Datos

Voy a usar una versión del mapa de México algo simplificada que pesa relativamente poco. Este shapefile puede descargarse aquí.

Tip: la versión más reciente de la cartografía de los municipios de México puede descargarse desde el Marco Geoestadístico Nacional de INEGI.

Leemos el archivo, asignamos la clave municipal de INEGI como el índice del GeoDataFrame y además reproyectamos las coordenadas usando el sistema de coordenadas Mexico ITRF92 / UTM zone 12N

mx = gpd.read_file('mapa_mexico/')\
        .set_index('CLAVE')\
        .to_crs(epsg=4485)
mx.head()
NOM_MUN NOMEDO CVE_EDO CVE_MUNI Area geometry
CLAVE
02004 Tijuana Baja California 02 004 1122.661145 POLYGON ((-73565.018 3602427.487, -73564.403 3...
02003 Tecate Baja California 02 003 3670.991923 POLYGON ((-38995.078 3617846.589, -31557.921 3...
02002 Mexicali Baja California 02 002 13119.275713 POLYGON ((48160.716 3621731.593, 58570.990 362...
02005 Playas de Rosarito Baja California 02 005 517.120801 POLYGON ((-70946.724 3594803.753, -70966.034 3...
26055 San Luis Rio Colorado Sonora 26 055 9033.770278 POLYGON ((127160.493 3587762.823, 127099.688 3...

Al GeoDataFrame anterior vamos a añadirle datos para que podamos graficar un mapa coroplético. Usaremos la API de DataMéxico para obtener los datos del Índice de Complejidad Económica (ECI). La propia definición del ECI de DataMexico es:

El Índice de Complejidad Económica (ECI) mide las capacidades productivas de una localidad (e.g. estado o municipalidad) a partir de la presencia de actividades (e.g. empleo, industrias o exportaciones) en esa y otras localidades. La complejidad económica de una localidad predice su nivel de ingreso, crecimiento económico, desigualdad, y emisiones de gases de efecto invernadero.

En general, un ECI más alto indica un mayor nivel de sofisticación en los productos que se producen en un municipio. Así se ven los datos:

url = 'https://api.datamexico.org/tesseract/cubes/complexity_eci/aggregate.jsonrecords?cuts%5B%5D=Latest.Latest.Latest.1&drilldowns%5B%5D=Geography+Municipality.Geography.Municipality&drilldowns%5B%5D=Date+Day.Date.Year&measures%5B%5D=ECI&parents=false&sparse=false'
data_eci = requests.get(url).json()
eci = pd.DataFrame(data_eci['data'])\
    .assign(CLAVE=lambda x: x['Municipality ID'].astype(str).str.zfill(5))\
    .set_index('CLAVE')

eci.head()
Municipality ID Municipality Year ECI
CLAVE
01001 1001 Aguascalientes 2020 2.976280
01002 1002 Asientos 2020 -0.430667
01003 1003 Calvillo 2020 0.420862
01004 1004 Cosío 2020 -0.464407
01005 1005 Jesús María 2020 3.049664

El índice está estandarizado con una media de 0 y varianza de 1 y tiene un rango entre -1.4 y 4.9.

eci['ECI'].describe()
count    2.142000e+03
mean     7.507756e-10
std      1.000000e+00
min     -1.399253e+00
25%     -6.253233e-01
50%     -2.969113e-01
75%      2.492735e-01
max      4.901706e+00
Name: ECI, dtype: float64

Revisamos su distribución:

eci['ECI'].hist(bins=100)
<AxesSubplot:>

Ahora anexamos los datos del ECI al GeoDataFrame:

mx = mx.join(eci, how='left')
mx.head()
NOM_MUN NOMEDO CVE_EDO CVE_MUNI Area geometry Municipality ID Municipality Year ECI
CLAVE
01001 Aguascalientes Aguascalientes 01 001 1168.762384 POLYGON ((1416489.577 2467700.472, 1417908.226... 1001.0 Aguascalientes 2020.0 2.976280
01002 Asientos Aguascalientes 01 002 547.762077 POLYGON ((1417043.958 2491681.240, 1417408.488... 1002.0 Asientos 2020.0 -0.430667
01003 Calvillo Aguascalientes 01 003 931.300088 POLYGON ((1347882.273 2454901.097, 1348002.307... 1003.0 Calvillo 2020.0 0.420862
01004 Cosio Aguascalientes 01 004 128.907513 POLYGON ((1397788.297 2509816.078, 1398009.089... 1004.0 Cosío 2020.0 -0.464407
01005 Jesus Maria Aguascalientes 01 005 499.207990 POLYGON ((1388272.165 2462097.533, 1389832.232... 1005.0 Jesús María 2020.0 3.049664

Creamos un par de GeoDataFrames que nos servirán más adelante:

  • oax contiene los datos solo para el estado de Oaxaca
  • edos es un GeoDataFrame con los polígonos de los estados del país.
oax = mx.query('CVE_EDO=="20"')
edos = mx.dissolve(by='CVE_EDO')

Y ahora visualicemos el ECI con un esquema de cortes por quintiles:

fig, ax = plt.subplots()
mx.plot(column='ECI', legend=True, scheme='quantiles', k=5, cmap='viridis_r', ax=ax)
fig.set_size_inches(12, 8)

Como se puede ver, algunas zonas del centro y sur del país son difíciles de distinguir por el gran número de municipios que hay.

Recuadro de acercamiento

Vamos a añadir un recuadro de acercamiento o detalle a la zona sur, señalando principalmente los municipios del estado de Oaxaca. Para eso usaremos la función inset_axes de Matplotlib. Esta función nos permite añadir un nuevo objeto axes a nuestra gráfica en cualquier parte del gráfico. La función recibe una lista con 4 valores que son: $[x, y, w, h]$:

  • $x$: posición en el eje x de la esquina inferior izquierda del recuadro
  • $y$: posición en el eje y de la esquina inferior izquierda del recuadro
  • $w$: ancho del recuadro
  • $h$: altura del recuadro

Todos estos valores vienen expresados como proporción del tamaño de la gráfica (proporción del ancho para el eje x, proporción del alto para el eje y). Veamos un ejemplo:

fig, ax = plt.subplots()
axins = ax.inset_axes([1.1, 0.5, 0.4, 0.4])
axins.set(xlim=(0.5, 0.7), ylim=(0.4, 0.6))
ax.indicate_inset_zoom(axins)
(<matplotlib.patches.Rectangle at 0x2f7a3bf22e0>,
 (<matplotlib.patches.ConnectionPatch at 0x2f7a3c269d0>,
  <matplotlib.patches.ConnectionPatch at 0x2f7a3c26d30>,
  <matplotlib.patches.ConnectionPatch at 0x2f7a3c26fa0>,
  <matplotlib.patches.ConnectionPatch at 0x2f7a0951550>))

Lo que hicimos fue añadir una figura (fig) y un objeto axes principal (ax). Luego añadimos el recuadro axins y limitamos el rango del eje $x$ e $y$ en los que se enfoca el recuadro. Por último, le indicamos al recuadro cuál es la gráfica principal sobre la que esta haciendo el zoom. Lo que falta es dibujar sobre los objetos ax y axins.

x = [0, 0.5, 0.6, 0.8]
y = [0, 0.4, 0.5, 0.3]
ax.plot(x, y, color='C0')
axins.plot(x, y, color='C0')
fig

Aunque es relativamente fácil añadir un recuadro a una gráfica, es un poco más difícil añadir el recuadro a un mapa. La dificultad adicional está en que en el resultado final los parámetros $x$, $y$, $w$ y $h$ no se comportan exactamente como sería esperado, sino que por algún motivo el recuadro queda en una posición un poco diferente. No podría asegurarlo, pero me parece que esto pasa porque al graficar un mapa Matplotlib modifica el tamaño de la gráfica y su relación de aspecto para ajustarlo a la información geográfica. La solución (muy poco satisfactoria) es jugar con varios valores de $x$, $y$, $w$ y $h$ hasta que la posición sea la adecuada.

En el siguiente ejemplo vamos a acercarnos al estado de Oaxaca con un recuadro usando los parámetros $x=0.82, y=0.05, w=0.9, h=0.9$ y poniendo como límite al eje x el rango (1820000, 2350000) y al eje y el rango (1780000, 2150000). Los valores (y su unidad de medida) de los ejes X e Y vienen dadas por el sistema de coordenadas geográficas que usamos.

fig, ax = plt.subplots()
# Añade recuadro
axins = ax.inset_axes([0.82, 0.05, 0.9, 0.9])

# Gráfica principal
mx.plot(column='ECI', legend=True, ax=ax, scheme='quantiles', k=5, cmap='viridis_r', legend_kwds={'loc': 'lower left'})
mx.boundary.plot(lw=0.05, color='k', ax=ax)

# Gráfica recuadro
mx.plot(column='ECI', legend=False, ax=axins, scheme='quantiles', k=5, cmap='viridis_r')
mx.boundary.plot(lw=0.25, color='k', ax=axins)
edos.boundary.plot(lw=1, color='red', ax=axins)

# limita área a mostrar
axins.set(ylabel='', xlabel='', xlim=(1820000, 2350000), ylim=(1780000, 2150000), xticks=[], yticks=[])

# Establece líneas del recuadro a la gráfica principal
ax.indicate_inset_zoom(axins)

# Tamaño de la gráfica final
fig.set_size_inches(12, 8)

El acercamiento resulta bastante bueno, ahora es posible identificar mejor los valores del ECI para el estado de Oaxaca y otros vecinos. Podemos añadir un segundo recuadro para ver los municipios del industrial estado de Guanajuato. El proceso es el mismo que el anterior, aunque ahora demorará más porque debe graficar más elementos.

fig, ax = plt.subplots()
# Añade recuadro 1
axins = ax.inset_axes([0.86, 0.05, 0.9, 0.9])
# Añade recuadro 2
axins2 = ax.inset_axes([0.2, -0.5, 0.5, 0.5])

# Gráfica principal
mx.plot(column='ECI', legend=True, ax=ax, scheme='quantiles', k=5, cmap='viridis_r', legend_kwds={'loc': 'lower left'})
mx.boundary.plot(lw=0.05, color='k', ax=ax)

# Gráfica recuadro 1
mx.plot(column='ECI', legend=False, ax=axins, scheme='quantiles', k=5, cmap='viridis_r')
mx.boundary.plot(lw=0.25, color='k', ax=axins)
edos.boundary.plot(lw=1, color='red', ax=axins)

# Gráfica recuadro 2
mx.plot(column='ECI', legend=False, ax=axins2, scheme='quantiles', k=5, cmap='viridis_r')
mx.boundary.plot(lw=0.25, color='k', ax=axins2)
edos.boundary.plot(lw=1, color='red', ax=axins2)

# limita área a mostrar recuadro 1 y 2
axins.set(ylabel='', xlabel='', xlim=(1820000, 2350000), ylim=(1780000, 2150000), xticks=[], yticks=[], title='Oaxaca')
axins2.set(ylabel='', xlabel='', xlim=(1425000, 1680000), ylim=(2229203, 2453767), xticks=[], yticks=[], title='Guanajuato')

# Elimina marco de la gráfica principal
ax.set_axis_off()
# Establece líneas de los recuadros a la gráfica principal
ax.indicate_inset_zoom(axins)
ax.indicate_inset_zoom(axins2)

# Tamaño de la gráfica final
fig.set_size_inches(12, 8)

Es notable la diferencia en la complejidad entre los municipios de Guanajuato y Oaxaca.

El pequeño recuadro al interior que señala el área del mapa principal que se está señalando no se nota mucho porque se confunde con el color de los polígonos. Para evitar eso, podemos recolorear el borde de los recuadros para que sea posible asociarlos usando el color. Esta idea la vi en un mapa creado por @ElenoAM 🤯 y me pareció incluso más limpia que tener varias líneas indicando el zoom.

Veamos ahora cómo queda tras quitar las líneas y ponerle colores al borde de los recuadros.

fig, ax = plt.subplots()
# Añade recuadro 1
axins = ax.inset_axes([0.86, 0.05, 0.9, 0.9])
# Añade recuadro 2
axins2 = ax.inset_axes([0.2, -0.5, 0.5, 0.5])

# Gráfica principal
mx.plot(column='ECI', legend=True, ax=ax, scheme='quantiles', k=5, cmap='viridis_r', legend_kwds={'loc': 'lower left'})
mx.boundary.plot(lw=0.05, color='k', ax=ax)

# Gráfica recuadro 1
mx.plot(column='ECI', legend=False, ax=axins, scheme='quantiles', k=5, cmap='viridis_r')
mx.boundary.plot(lw=0.25, color='k', ax=axins)
edos.boundary.plot(lw=1, color='red', ax=axins)

# Gráfica recuadro 2
mx.plot(column='ECI', legend=False, ax=axins2, scheme='quantiles', k=5, cmap='viridis_r')
mx.boundary.plot(lw=0.25, color='k', ax=axins2)
edos.boundary.plot(lw=1, color='red', ax=axins2)

# limita área a mostrar recuadro 1 y 2
axins.set(ylabel='', xlabel='', xlim=(1820000, 2350000), ylim=(1780000, 2150000), xticks=[], yticks=[], title='Oaxaca')
axins2.set(ylabel='', xlabel='', xlim=(1425000, 1680000), ylim=(2229203, 2453767), xticks=[], yticks=[], title='Guanajuato')

# Elimina marco de la gráfica principal
ax.set_axis_off()
# Establece líneas de los recuadros a la gráfica principal
ax.indicate_inset_zoom(axins)
ax.indicate_inset_zoom(axins2)

# Parametros para cambiar colores
color_insets = ['red', 'blue']
color_index = 0

# colorea los cuadros interiores
for p in ax.patches:
    if isinstance(p, plt.matplotlib.patches.Rectangle):
        p.set_edgecolor(color_insets[color_index])
        p.set_linewidth(2)
        color_index += 1
    else:
        # Esconde las líneas hacia la gráfica principal
        p.set_visible(False)

# Colorea el marco del recuadro 1 y aumenta grosor de línea
for sp in axins.spines.values():
    sp.set_color(color_insets[0])
    sp.set_linewidth(2)

# Colorea el marco del recuadro 2 y aumenta grosor de línea
for sp in axins2.spines.values():
    sp.set_color(color_insets[1])
    sp.set_linewidth(2)
    
# Tamaño final
fig.set_size_inches(12, 8)

Ahora es más fácil identificar el área que se está acercando.

Recuadro de contexto

En un recuadro de contexto el mapa principal corresponde a la región que vamos a analizar con detalle y usamos el recuadro para resaltar la posición de esa región con respecto al país. Para lograr este efecto no usamos la función ax.inset_axes(), sino una muy parecida que es fig.add_axes() y recibe los mismos argumentos $x, y, w, h$. La diferencia es que ahora no tendremos las líneas desde el recuadro hasta el mapa principal.

fig, ax = plt.subplots()
# Mapa principal
oax.plot(column='ECI', legend=True, scheme='quantiles', k=5, cmap='viridis_r', ax=ax)
# Borde de los municipios
oax.boundary.plot(linewidth=0.5, color='gray', ax=ax)
# Inserta recuadro
ax_mex = fig.add_axes([0.91, 0.11, 0.35, 0.35], )
# Dibuja los estados del país en el recuadro
edos.boundary.plot(color='gray', ax=ax_mex)
# Resalta Oaxaca con el color rojo
edos.query('CVE_EDO=="20"').plot(color='red', ax=ax_mex)
# Establece tamaño final
fig.set_size_inches(12, 8)

Escalas

Cuando graficamos mapas con distintos niveles de acercamiento/alejamiento es importante acompañar cada mapa con una barra de escala que nos permita saber cómo comparar longitudes. Para Matplotlib existe la excelente librería matplotlib-scalebar para agregar estas escalas y afortunadamente funciona bien con mapas de Geopandas. A continuación añadimos las escalas para cada mapa.

ax.set(xticks=[], yticks=[], title='Oaxaca')
ax_mex.set(xticks=[], yticks=[], title='México')
# Barras de escala
scalebar_oax = ScaleBar(1, "m", length_fraction=0.2, location='lower left', )
ax.add_artist(scalebar_oax)
scalebar_mex = ScaleBar(1, "m", length_fraction=0.25, location='lower left')
ax_mex.add_artist(scalebar_mex)
fig

De esta forma podemos reconocer mejor las dimensiones de Oaxaca y sus municipios con respecto a todo el país.

Recuadro de áreas no contiguas

Este tipo de recuadro sirve para mostrar juntas áreas que no están cerca o es muy difícil ver juntas en su ubicación normal. El primer caso que viene a mi cabeza es el del departamento de San Andrés y Providencia en Colombia 🟨🟨🟦🟥. San Andrés y Providencia son unas diminutas islas ubicadas a 774 km del noreste de Colombia y que conforman uno de los 32 departamentos del país. En los mapas de división política-administrativa del país estas islas siempre aparecen como recuadros al lado de la parte continental, con una escala mucho mayor para poder ser visibles.

Para visualizarlas usemos este mapa de Colombia con la división por departamentos que también puede descargarse desde el Marco Geostadístico Nacional. Usaremos el sistema de coordenadas MAGNA-SIRGAS / Colombia Bogota zone.

col = gpd.read_file('mapa_colombia/').to_crs(epsg=3116)
col.head()
DPTO_CCDGO DPTO_CNMBR DPTO_NANO_ DPTO_CACTO DPTO_NANO SHAPE_AREA SHAPE_LEN geometry
0 18 CAQUETÁ 1981 Ley 78 del 29 de Diciembre de 1981 2020 7.318485 21.384287 POLYGON ((909199.769 818940.314, 909214.572 81...
1 19 CAUCA 1857 15 de junio de 1857 2020 2.534419 13.950263 POLYGON ((735237.006 860162.569, 735286.223 86...
2 86 PUTUMAYO 1991 Articulo 309 Constitucion Politica de 1991 2020 2.107965 12.707922 POLYGON ((711344.106 654182.557, 711400.239 65...
3 76 VALLE DEL CAUCA 1910 Decreto No 340 de 16 de Abril de 1910 2020 1.679487 12.650870 MULTIPOLYGON (((777973.333 1049936.014, 777964...
4 94 GUAINÍA 1991 Articulo 309 Constitucion Politica de 1991 2020 5.747937 21.179051 POLYGON ((1712400.898 927096.060, 1712774.632 ...

Al graficar el mapa del país a duras penas podremos notar las islas en la parte superior derecha.

fig, ax = plt.subplots()
col.plot(ax=ax)
fig.set_size_inches(10, 10)

Para hacer un típico mapa político administrativo, vamos a añadir dos objetos axes con fig.add_axes() en los que graficaremos a cada isla. Para encontrar las coordenadas de las islas lo mejor es usar el mapa de http://epsg.io/map#srs=3116&x&y&z=10&layer=streets. También agregamos las barras de escala para dimensionar mejor lo pequeñas que son.

fig, ax = plt.subplots()
# Añade ax para San Andrés
ax_san = fig.add_axes([0.25, 0.75, 0.05, 0.1])
# Añade ax para Providencia
ax_prov = fig.add_axes([0.29, 0.77, 0.08, 0.08])

# Grafica Parte continental
col.plot(ax=ax)
# Grafica San Andrés
col.plot(ax=ax_san)
# Grafica Providencia
col.plot(ax=ax_prov)

# Limita los valores del eje x e y para cada gráfica
ax.set(ylim=(0, 1_900_000), xlim=(400_000, 1_850_000), yticks=[], xticks=[])
ax_san.set(ylim=(1_883_000, 1_902_000), xlim=(163_000, 172_000), yticks=[], xticks=[])
ax_prov.set(ylim=(1_975_000, 1_990_000), xlim=(204_044, 211_100), yticks=[], xticks=[])

# Títulos de las gráficas
ax.set_title('Colombia', fontsize=20)
ax_san.set_title('San\n Andrés', fontsize=7)
ax_prov.set_title('Providencia', fontsize=7)

# Añade barras de escala
scalebar_col = ScaleBar(1, "m", length_fraction=0.2, location='lower left', )
ax.add_artist(scalebar_col)
scalebar_san = ScaleBar(1, "m", length_fraction=1, location='upper center', font_properties={'size': 7})
ax_san.add_artist(scalebar_san)
scalebar_prov = ScaleBar(1, "m", length_fraction=0.8, location='upper center',  box_alpha=0, font_properties={'size': 7})
ax_prov.add_artist(scalebar_prov)

# Tamaño final de la gráfica
fig.set_size_inches(10, 10)