Introducción

Los mapas coropléticos son una de las formas más comunes de presentar información geográfica. En QGIS se pueden hacer mapas muy sencillos o también se puede lograr diseños más complejos, sin embargo, una desventaja de usar este software para crear un mapa coroplético es que normalmente es un proceso lento que requiere de muchos pasos (comparado con hacerlo en R o Python, por ejemplo). Si esta tarea se realiza manualmente entonces se requiere de mucho esfuerzo para hacer varios mapas en serie. Y si ya hiciste los mapas pero solo quieres cambiar un elemento, como por ejemplo, la paleta de colores, entonces habría que editar nuevamente cada uno de ellos.

Afortunadamente con PyQGIS es posible crear mapas usando código de Python, lo que puede reducir la carga de trabajo y automatizar ciertos pasos que nos ahorren tiempo. Aun así desde ahora advierto: hacer mapas coropléticos, y en general, manejar la simbología de capas con PyQGIS no es una tarea fácil. Se requiere de un buen conocimiento sobre cómo funciona PyQGIS. Y aún sabiendo de PyQGIS, hacer un mapa coroplético conlleva muchos pasos, por lo que hay que evaluar bien si vale la pena invertir tiempo en programar todo un proceso para generar estos mapas. Yo por ejemplo, lo encontré útil porque estaba haciendo un proceso iterativo en el que editaba manualmente algunos shapefiles y luego quería ver rápidamente el resultado en un mapa coroplético. Como necesitaba hacer esto muchas veces, valió la pena hacer la inversión de tiempo en automatizar los mapas, pero si en tu caso solo necesitas hacer pocos mapas, yo recomendaría hacerlo manual o usar otras alternativas más rápidas.

Esta entrada será larga, porque se requiere entender muchos conceptos para poder llegar a hacer algo tan sencillo como un mapa coroplético en QGIS. Y ya aprovechando la explicación, voy a exponer en general el tema de la simbología de capas y a mostrar varios ejemplos. Esta entrada de hecho la hago más para mí cuando vuelva a necesitar esto en el futuro.

Si no conoces de PyQGIS, puedes revisar la entrada que escribí al respecto: Introducción a PyQGIS (Python + QGIS)

Configuración inicial

Esta entrada la escribo desde un ambiente de conda (en Windows 10), tal como expliqué en la entrada Integración de PyQGIS con Jupyter Lab usando Anaconda, aunque el código debería funcionar igual en la consola de Python integrada dentro de QGIS (si no tiene mucha esperiencia, recomiendo usar la consola integrada). Si también quieres correr este código desde un ambiente de conda, entonces debes correr el siguiente código al inicio de tu script:

from qgis.core import QgsApplication
from IPython.display import Image
qgs = QgsApplication([], False)
qgs.initQgis()

Warning: si usas la consola de Python dentro de QGIS entonces no es necesario ejecutar el código anterior

Con el siguiente código vamos a importar las librerías que necesitamos y además voy a definir algunas funciones auxiliares que me ayudarán a exportar los resultados y poder incluirlos en este notebook. Más adelante explicaremos lo que hace cada función, por ahora simplemente basta con definirlas.

import os
from qgis.core import *
from PyQt5.QtGui import QColor
from PyQt5.QtCore import QSize, Qt
from qgis.PyQt.QtGui import QFont

dir_data = "datos"


def create_layout_with_map(layer, scale_factor=1):

    layout = QgsPrintLayout(proyecto)
    layout.initializeDefaults()

    page = layout.pageCollection().page(0)
    pagesize = page.pageSize()

    mapa = QgsLayoutItemMap(layout)
    extent = layer.extent()
    # Set map item position and size (by default, it is a 0 width/0 height item placed at 0,0)
    margin_y = 10
    layout_height = pagesize.height() - 2 * margin_y
    layout_width = min(layout_height * (extent.width() / extent.height()), 2 * pagesize.width())
    margin_x = (pagesize.width() / 2) - (layout_width / 2)
    mapa.attemptMove(QgsLayoutPoint(margin_x, margin_y, QgsUnitTypes.LayoutMillimeters))
    mapa.attemptResize(QgsLayoutSize(layout_width, layout_height, QgsUnitTypes.LayoutMillimeters))
    mapa.setExtent(extent)
    mapa.setScale(mapa.scale() * scale_factor)
    layout.addLayoutItem(mapa)
    return layout


def add_legend_to_layout(layout, title='', pos_x=15, pos_y=10):
    mapa = layout.referenceMap()
    legend = QgsLayoutItemLegend(layout)
    legend.setLinkedMap(mapa)  # map is an instance of QgsLayoutItemMap
    layout.addLayoutItem(legend)
    legend.attemptMove(QgsLayoutPoint(pos_x, pos_y, QgsUnitTypes.LayoutMillimeters))
    legend.setTitle(title)


def export_image(layout, img_path: str, dpi: int = 100):
    exporter = QgsLayoutExporter(layout)
    exp_sett = QgsLayoutExporter.ImageExportSettings()
    exp_sett.dpi = dpi
    exporter.exportToImage(img_path, exp_sett)

Tip: si trabajas desde un Jupyter Notebook, puedes también definir la siguiente función y de esta manera puedes revisar cómo queda el resultado después de ejecutar alguna instrucción con show_image(layer), donde layer es la capa que queremos visualizar.

def show_image(layer, scale_factor=1, dpi=50, legend_title=''):
    os.makedirs('images', exist_ok=True)
    img_path = 'images/simbologia.png'
    layout = create_layout_with_map(layer, scale_factor=scale_factor)
    add_legend_to_layout(layout, title=legend_title)
    export_image(layout, img_path, dpi=dpi)
    return Image(img_path, width=600)

Vamos cargar una capa que contiene las colonias de la CDMX. Los datos que uso se pueden descargar desde aquí. Tras ejecutar el código deberíamos observar simplemente la capa ya cargada. Todos los polígonos se verán de un mismo color.

proyecto = QgsProject.instance()
proyecto.removeAllMapLayers()

nombre_layer_colonias = "colonias_cdmx"
layer_colonias = QgsVectorLayer(f"{dir_data}/colonias_iecm_2019/mgpc_2019.shp", nombre_layer_colonias)
proyecto.addMapLayer(layer_colonias)
<QgsVectorLayer: 'colonias_cdmx' (ogr)>

Resumen mapas coropléticos

Esta entrada es muy larga y requiere tiempo procesar todo lo que aquí trato de explicar. Si lo que buscas concretamente es cómo hacer un mapa coroplético a partir de una capa que contiene una variable contínua y además añadir la leyenda correspondiente y que se exporte a una imagen en formato .png, entonces voy a resumir esa tarea en la siguiente función (además de las funciones definidas anteriormente):

def create_choropleth(layer: QgsVectorLayer, value_field: str, num_classes: int = 5, classification_method='Quantile',
                      ramp_name: str = 'Spectral', img_path: str = 'image.png'):
    """
    Genera un mapa coroplético de la capa `layer`, usando la variable `value_field`. El resultado es exportado a `img_path`
    `layer`: es la capa que queremos representar en el mapa.
    `value_field`: es el nombre la variable continua que nos interesa representar en el mapa. Los polígonos tendrán un color asignado de acuerdo con el valor que tengan en esta variable.
    `num_classes`: es el número de categorías en que queremos dividir los valores de la variable en value_field.
    `classification_method`: es el método con el que se formarán los intervalos para conformar las categorías señaladas por `num_classes`. Las opciones pueden ser: 'Quantile', 'EqualInterval', 'Logarithmic', 'Jenks', 'Pretty', 'StdDev'
    `ramp_name`: es el nombre de la paleta de colores que se usará.
    `img_path`: ruta y nombre del archivo donde se guardará la imagen generada.
    """

    proyecto.removeAllMapLayers()
    proyecto.addMapLayer(layer)
    default_style = QgsStyle().defaultStyle()
    color_ramp = default_style.colorRamp(ramp_name)

    renderer = QgsGraduatedSymbolRenderer()
    renderer.setClassAttribute(value_field)
    classification_method_instance = QgsClassificationMethodRegistry().method(classification_method)
    renderer.setClassificationMethod(classification_method_instance)
    renderer.updateClasses(layer, num_classes)
    renderer.updateColorRamp(color_ramp)
    layer.setRenderer(renderer)

    layout = create_layout_with_map(layer)
    add_legend_to_layout(layout, title=value_field)
    export_image(layout, img_path)

La función anterior tiene 6 parámetros que vienen explicados en la descripción de la misma. Para nuestro ejemplo, vamos a representar la variable POB2010 que es la población de cada colonia. Escogimos representar en 5 clases, usando el método de clasificación de Jenks y la paleta de colores Spectral. Así se ejecuta y se ve el resultado:

layer_colonias = QgsVectorLayer(f"{dir_data}/colonias_iecm_2019/mgpc_2019.shp", nombre_layer_colonias)
img_path = 'images/choropleth.png'

create_choropleth(layer=layer_colonias, value_field='POB2010', num_classes=5, classification_method='Jenks',
                  ramp_name='Spectral', img_path=img_path)
  • Los campos disponibles de la capa layer_colonias son:
layer_colonias.fields().toList()
[<QgsField: ENT (String)>,
 <QgsField: CVEDT (String)>,
 <QgsField: NOMDT (String)>,
 <QgsField: DTTOLOC (String)>,
 <QgsField: CVEUT (String)>,
 <QgsField: NOMUT (String)>,
 <QgsField: POB2010 (String)>,
 <QgsField: ID (Integer)>]
  • Las opciones disponibles para ramp_name se puede consultar con el siguiente comando:
QgsStyle().defaultStyle().colorRampNames()
['Blues',
 'BrBG',
 'BuGn',
 'BuPu',
 'GnBu',
 'Greens',
 'Greys',
 'Inferno',
 'Magma',
 'OrRd',
 'Oranges',
 'PRGn',
 'PiYG',
 'Plasma',
 'PuBu',
 'PuBuGn',
 'PuOr',
 'PuRd',
 'Purples',
 'RdBu',
 'RdGy',
 'RdPu',
 'RdYlBu',
 'RdYlGn',
 'Reds',
 'Spectral',
 'Viridis',
 'YlGn',
 'YlGnBu',
 'YlOrBr',
 'YlOrRd']
  • Las opciones para classification_method son:
list(QgsClassificationMethodRegistry().methodNames().values())
['Quantile', 'EqualInterval', 'Logarithmic', 'Jenks', 'Pretty', 'StdDev']

Usemos la función en otro ejemplo con otros inputs:

layer_colonias = QgsVectorLayer(f"{dir_data}/colonias_iecm_2019/mgpc_2019.shp", nombre_layer_colonias)
create_choropleth(layer=layer_colonias, value_field='ID', num_classes=10, classification_method='Pretty',
                  ramp_name='Plasma', img_path=img_path)