Introducción

Últimamente he estado trabajando mucho con QGIS, el software SIG de código libre indispensable para cualquier analista de datos geográficos. El trabajo que realizaba era muy manual (crear y editar polígonos), pero había ciertos pasos que sentía que podía automatizar, así que decidí aprender al menos un poco de PyQGIS para hacer las tareas repetitivas.

PyQGIS es la librería de Python para ejecutar tareas dentro de QGIS. Si has usado QGIS seguro que en más de una vez has visto el logo de Python en la barra de herramientas. Python y QGIS tienen una gran integración y aprender a usar PyQGIS es una inversión de tiempo con rendimientos garantizados. Pese a ello, me parece que no es muy fácil iniciar con PyQGIS, incluso si uno ya sabe Python. PyQGIS y su documentación parece que están diseñados más hacia desarrolladores de plugins que hacia usuarios finales. Si uno viene de utilizar herramientas de análisis geográfico más sencillas como GeoPandas o SF seguramente PyQGIS parece muy burocrático o complicado en un principio, pero como en todo, se trata de entender su lógica y poco a poco uno va descubriendo que también es una herramienta muy completa y muy práctica.

En esta entrada me gustaría mostrar la lógica básica para que alguien que ya tiene conocimentos de la sintaxis de Python y ha trabajado con QGIS pueda empezar a automatizar procesos sencillos con PyQGIS. En particular voy a mostrar un par de ejemplos prácticos aplicando algoritnmos de validación y procesamiento. Aclaro que no soy un experto en QGIS, pero quiero compartir mi aprendizaje al iniciar porque sé que puede haber otras personas en la misma situación.

La consola de Python en QGIS

Voy a trabajar desde la interfaz gráfica de QGIS, por lo que obviamente hay que tener instalado QGIS en la computadora. En otra entrada mostraré cómo instalar QGIS desde Anaconda para trabajar usando JupyterLab.

Una vez instalado y descargado QGIS, abrimos el programa e iniciamos un nuevo proyecto. En la barra de herramientas encontraremos el símbolo de Python. Al hacer clic se abrirá una consola de Python dentro de la ventana de QGIS. Si no ves el ícono, puedes acceder usando el menú en la pestaña complementos >> Consola de Python

Como se ve en la siguiente imagen, la consola se divide en dos secciones, una es la ventana de resultados y la otro es la terminal de Python.

En la terminal puedes escribir código de Python directamente y este se mostrará en la ventana de resultados, por ejemplo, si escribimos en la terminal:

print("¡Hola mundo!")

en la ventana de resultados se imprimirá el comando que escribimos y su resultado

Escrbir código en la terminal no es muy práctico, por eso es mejor usar el editor de scripts, al que se puede acceder al hacer clic en el ícono superior de una hoja con un lápiz. También (me parece) es mejor trabajar en una ventana separada, para ello hacemos clic en el ícono de la parte superior derecha con dos cuadritos (al lado de la X).

En la nueva configuración, el editor de scripts aparece a la derecha y la ventana de resultados, junto con la terminal, del lado izquierdo. Podemos escribir código en el editor y ejecutarlo haciendo clic en el ícono superior en forma de triángulo verde (o con el atajo de teclado ctrl + shift + e), el resultado se imprime en la ventana de resultados.

Se puede guardar el script con el ícono de disco, o se puede crear uno nuevo al hacer clic en el símbolo $+$

Básicos de PyQGIS

Dentro de la consola de Python tenemos acceso a las librerías con las funciones importantes para poder manipular elementos de QGIS. Vamos a crear un nuevo script que se llame intro_pyqgis_01.py. En el importamos las funciones (constructores, para ser más precisos) QgsProject, QgsVectorLayer del módulo qgis.core.

from qgis.core import QgsProject, QgsVectorLayer

Como sabemos, en QGIS se usan proyectos que contienen capas (layers). Las capas normalmente hacen referencia a archivos con datos geoespaciales, típicamente Shapefiles o GeoJSON, aunque también puede haber capas de otros tipos. Bueno, esa misma estructura la debemos tener en cuenta en PyQGIS. Debemos trabajar sobre un objeto tipo proyecto que contiene objetos de tipo capas.

Vamos a iniciar manualmente un nuevo proyecto de QGIS (proyecto >> Nuevo) y lo guardamos en nuestro equipo con el nombre de intro_pyqgis.

Ahora, vamos al editor de script y creamos una instancia de proyecto a partir del proyecto que actualmente está activo en QGIS:

from qgis.core import QgsProject, QgsVectorLayer
proyecto = QgsProject.instance()

Al elecutar lo anterior no veremos que nada cambie, solo definimos el objeto proyecto. Podemos ecribir la siguiente línea en la terminal para verificar el archivo de proyecto que estamos usando:

print(proyecto.fileName())

Que en mi caso arroja:

> >> 'C:/Users/juan.santos/OneDrive/intro_pyqgis/intro_pyqgis.qgz'

Ahora hay que agregarle datos (capas). Para eso usamos la función QgsVectorLayer, que recibe dos argumentos: ruta del archivo geográfico y nombre de la capa.

Los datos que usaremos son publicados por el portal de datos abiertos de la Ciudad de México (CDMX):

Como los datos pueden cambiar en el futuro, dejo aquí para descarga los que yo usé en particular para este post.

Ahora sí vamos a añadir una capa con los datos de las colonias de la CDMX que se va a llamar "colonias_cdmx".

## intro_pyqgis_01.py
from qgis.core import QgsProject, QgsVectorLayer

# Ruta a la carpeta datos
dir_data = "datos"

proyecto = QgsProject.instance()

# Nombre de la capa
nombre_layer_colonias = "colonias_cdmx"

# crea la capa
layer_colonias = QgsVectorLayer(f"{dir_data}/colonias_iecm_2019/mgpc_2019.shp", nombre_layer_colonias)

# Añade la capa
proyecto.addMapLayer(layer_colonias)

Ahora sí veremos que en la ventana principal de QGIS se carga el mapa de colonias.

¡Y así de fácil es cargar datos! Sin embargo, debemos de tener cuidado porque si corremos el mismo script varias veces se cargará la misma capa muchas veces, ya que QGIS permite cargar múltiples veces el mismo archivo.

Para evitar eso, vamos a modificar el script e incluir una pequeña condición que primero verifica si ya existe en el proyecto una capa con el mismo nombre. En caso de que no esxista, entonces carga la capa, si ya existe, entonces no hace nada más. De esta forma podemos ejecutar el archivo intro_pyqgis_01.py tantas veces queramos sin cargar la misma capa.

## intro_pyqgis_01.py
from qgis.core import QgsProject, QgsVectorLayer

# Ruta a la carpeta datos
dir_data = "datos"

proyecto = QgsProject.instance()

# Nombre de la capa
nombre_layer_colonias = "colonias_cdmx"

# crea la capa
layer_colonias = QgsVectorLayer(f"{dir_data}/colonias_iecm_2019/mgpc_2019.shp", nombre_layer_colonias)

# verifica si la capa ya existe
if not proyecto.mapLayersByName(nombre_layer_colonias):
    # Si no existe, Añade la capa
    proyecto.addMapLayer(layer_colonias)

Pues estos (QgsProject y QgsVectorLayer) me parece que son los tipos de objeto más básicos que hay que entender para empezar en PyQGIS (y saber que hay muchos, muchos más tipos de objetos).

Como todo objeto en Python, los proyectos y capas tienen métodos y atributos. De hecho, tienen muchos métodos y atributos, por lo que es dificil saber qué es lo que uno necesita. Ya vimos por ejemplo, un método para obtener el nombre del archivo del proyecto que estamos usando, y también cómo añadir una capa a un proyecto. Si queremos saber todos los métodos disponibles para un objeto projecto podemos consultar la documentación.

De ahí podemos obtener algunos datos útiles para trabajar:

print("Número de capas en el proyecto", proyecto.count())
print("CRS del proyecto", proyecto.crs().description())
print("Directorio del proyecto", proyecto.homePath())

resultado:

Número de capas en el proyecto 1
CRS del proyecto WGS 84 / UTM zone 14N
Directorio del proyecto C:/Users/juan.santos/OneDrive/intro_pyqgis

Igualmente, de las capas podemos encontrar cosas muy interesantes, como el número de filas y columnas que tienen los datos asociados, o el sistema de coordenadas. Consulta aquí la documentación de los objetos QgsVectorLayer

print("Número de filas:", layer_colonias.featureCount())
print("Número de columnas:", len(layer_colonias.fields()))
print("Sistema de coordenadas:", layer_colonias.crs().description())

resultado:

Número de filas: 1815
Número de columnas: 8
Sistema de coordenadas: WGS 84 / UTM zone 14N

Así que ya sabemos que hay 1815 colonias incluídas en este archivo y que tiene 8 columnas. ¿Cuál es el nombre de las columnas y de qué tipo son?:

for f in layer_colonias.fields():
    print(f.name(), f.typeName())

resultado:

ENT String
CVEDT String
NOMDT String
DTTOLOC String
CVEUT String
NOMUT String
POB2010 String
ID Integer

También podemos obtener los datos de una fila en particular de la capa con el método .getFeature(). El resultado es un objeto tipo QgsFeature. En este caso por ejemplo obtenemos la fila en la posición 11 del archivo y vemos los datos que contiene.

feature = layer_colonias.getFeature(10)
print(feature)
print(feature.attributes())

resultado:

<qgis._core.QgsFeature object at 0x000001A94F2EBB80>
['9', '2', 'AZCAPOTZALCO', '05', '02-013', 'CUITLÃ\x81HUAC 1 Y 2 (U HAB)', '2449', 11]

Y así, podemos seguir obteniendo más objetos. Por ejemplo del feature anterior podemos obtener la geometría, que es un objeto tipo QgsGeometry: Multipolygon. De la geometría, podemos obtener el centroide, que es un QgsGeometry: Point

geom = feature.geometry()
print(geom.centroid())

resultado:

<QgsGeometry: Point (482445.29628468106966466 2153086.92174857435747981)>

y así le vamos entendiendo de a poco, hasta parece sencillo.

Es muy interesante porque conociendo bien esta lógica podemos integrar muy bien QGIS con las demás herramientas de Python. Por ejemplo, solo por diversión, vamos a hacer una gráfica de la colonia que seleccionamos usando matplotlib, que ya viene instalado dentro del entorno de QGIS.

import matplotlib.pyplot as plt
import json

coords = json.loads(geom.asJson())["coordinates"][0][0]
coords_centroide = json.loads(geom.centroid().asJson())["coordinates"]
x, y = [[xy[n] for xy in coords] for n in [0, 1]]
ax = plt.subplot()
ax.plot(x, y)
ax.plot([coords_centroide[0]], [coords_centroide[1]], 'bo')
ax.text(x=coords_centroide[0], y=coords_centroide[1], s=feature.attribute("NOMUT"))
ax.figure.show()

Así queda el script ìntro_pyqgis_01.py completo:

# ìntro_pyqgis_01.py
from qgis.core import QgsProject, QgsVectorLayer

proyecto = QgsProject.instance()

dir_data = "datos"
nombre_layer_colonias = "colonias_cdmx"
layer_colonias = QgsVectorLayer(f"{dir_data}/colonias_iecm_2019/mgpc_2019.shp", nombre_layer_colonias)
if not proyecto.mapLayersByName(nombre_layer_colonias):
    proyecto.addMapLayer(layer_colonias)

print("Número de capas en el proyecto:", proyecto.count())
print("CRS del proyecto:", proyecto.crs().description())
print("Directorio del proyecto:", proyecto.homePath())

print("Número de filas:", layer_colonias.featureCount())
print("Número de columnas:", len(layer_colonias.fields()))
print("Sistema de coordenadas:", layer_colonias.crs().description())

for f in layer_colonias.fields():
    print(f.name(), f.typeName())
#
feature = layer_colonias.getFeature(10)
print(feature)
print(feature.attributes())
#
geom = feature.geometry()
print(geom.centroid())

import matplotlib.pyplot as plt
import json
coords = json.loads(geom.asJson())["coordinates"][0][0]
coords_centroide = json.loads(geom.centroid().asJson())["coordinates"]
x, y = [[xy[n] for xy in coords] for n in [0, 1]]
ax = plt.subplot()
ax.plot(x, y)
ax.plot([coords_centroide[0]], [coords_centroide[1]], 'bo')
ax.text(x=coords_centroide[0], y=coords_centroide[1], s=feature.attribute("NOMUT"))
ax.figure.show()

Operaciones con/entre capas

Uno de los usos que le podemos dar a PyQGIS, y en particular fue la razón por la que lo empecé a usar, es para automatizar operaciones con o entre capas.

Por ejemplo, yo estaba creando polígonos manualmente y cada tanto tenía que revisar que no tuvieran errores. Inicialmente hacía esto siguiendo la ruta del menú Vectorial >> Herramientas de Geometría >> Comprobar Validez y seleccionaba las opciones para que el algoritmo se ejecutara. Todo esto era repetitivo y además era medio molesto, porque como resultado de la ejecución siempre añadía 3 capas al proyecto, que me tocaba retirar manualmente. Entonces fue cuando descubrí que podía programar la operación de validación y así solo tenía que ejecutar el script cuando lo necesitara. Luego seguí automatizando otras operaciones que eran repetitivas.

Casi toda operación con capas que se puede hacer usando las opciones del menú de QGIS se puede hacer también usando PyQGIS de una manera relativamente sencilla. Podemos ver todas las operaciones/algoritmos disponibles en el menú Procesos > Caja de herramientas. En la caja de búsqueda podemos buscar alguna operación, por ejemplo, si lo que me interesa es calcular los centroides de una capa escribo centroide y obtengo los algoritmos relacionados:

Al pasar el cursos por encima de cualquier algoritmo aparece una pequeña nota con el ID del algoritmo en cuestión (native:centroids), lo que nos permitirá identificarlo para trabajar con PyQGIS.

Centroides

Para usar estos algoritmos necesitamos importar la librería processing. Con el ID del algoritmo podemos consultar en la terminal la información sobre los parámetros que espera recibir. Creamos un nuevo script con el nombre intro_pyqgis_02.py y ponemos las siguientes líneas:

# intro_pyqgis_02.py
import processing
processing.algorithmHelp("native:centroids")

que imprimen lo siguiente:

Centroides (native:centroids)

Este algoritmo crea una capa de puntos nueva, con puntos que representan el centroide de las geometrías de la capa de entrada.

Los atributos asociados a cada punto de la capa de salida son los mismos asociados a los objetos originales.


----------------
Input parameters
----------------

INPUT: Capa de entrada

    Parameter type: QgsProcessingParameterFeatureSource

    Accepted data types:
        - str: ID de la capa
        - str: Nombre de capa
        - str: Fuente de la capa
        - QgsProcessingFeatureSourceDefinition
        - QgsProperty
        - QgsVectorLayer

ALL_PARTS: Crear centroide para cada parte

    Parameter type: QgsProcessingParameterBoolean

    Accepted data types:
        - bool
        - int
        - str
        - QgsProperty

OUTPUT: Centroides

    Parameter type: QgsProcessingParameterFeatureSink

    Accepted data types:
        - str: archivo vectorial de destino, ej. 'd:/test.shp'
        - str: 'memory:' para guardar el resultado en una capa temporal en memoria
        - str: utilizando prefijo ID del proveedor vectorial y URI de destino, por ejemplo 'postgres:...' para guardar el resultado en una tabla PostGIS
        - QgsProcessingOutputLayerDefinition
        - QgsProperty

----------------
Outputs
----------------

OUTPUT:  <QgsProcessingOutputVectorLayer>
    Centroides

Con esta información sabemos que hay tres parámetros que puede recibir como input este algoritmo:

  • INPUT: Es la capa de la que se van a calcular los centroides.
  • ALL_PARTS: Indica si se va a crear centroide para cada parte, en caso de que se trate de multipolígonos
  • OUTPUT: dónde va a guardar el resultado. Puede ser un archivo o una capa temporal.

Algunos parámetros pueden ser opcionales. Estos se van a pasar a PyQGIS en un diccionario.

Así es como podemos ejecutar este algoritmo sobre la capa de colonias_cdmx:

# intro_pyqgis_02.py
from qgis.core import QgsProject, QgsVectorLayer
import processing

proyecto = QgsProject.instance()

dir_data = "datos"
nombre_layer_colonias = "colonias_cdmx"
if not proyecto.mapLayersByName(nombre_layer_colonias):
    layer_colonias = QgsVectorLayer(f"{dir_data}/colonias_iecm_2019/mgpc_2019.shp", nombre_layer_colonias)
    proyecto.addMapLayer(layer_colonias)

params = {
    'INPUT' : nombre_layer_colonias,
    'ALL_PARTS' : False,
    'OUTPUT' :'TEMPORARY_OUTPUT' 
}

out = processing.run("native:centroids", params)

Que desafortunadamente nos arroja el siguiente error:

_core.QgsProcessingException: Feature (11) from “colonias_cdmx” has invalid geometry. Please fix the geometry or change the Processing setting to the “Ignore invalid input features” option.

Esto debido a que la capa de colonias tiene errores en su creación :(. En un caso donde la capa no tuviera errores el algoritmo corriera sin problema y el resultado se guadaría en out.

Ahora usemos el algoritmo qgis:checkvalidity para verificar qué está mal con esta capa.

Validación

No voy a ahondar mucho en los detalles de los INPUTS de este algoritmo, pero pueden revisarlo con

processing.algorithmHelp("qgis:checkvalidity")

Solo notemos que estoy guardando el resultado en la variable out_validacion y que verifico si el algoritmo detecta características inválidas en la geometría. Si es así, entonces verifico si la capa resultante ya está en el proyecto. Si es así, entonces remuevo la capa existente. Luego, añado la capa que contiene los errores al proyecto. Esta lógica de remover y volver a añadir la capa se justifica porque en teoría yo esperaría que si detecto errores tendría que editar la capa original y luego ir verificando si el error se solucionó. Creamos un nuevo script con el nombre intro_pyqgis_03.py y ponemos las siguientes líneas:

# intro_pyqgis_03.py
from qgis.core import QgsProject, QgsVectorLayer
import processing

proyecto = QgsProject.instance()

dir_data = "datos"
nombre_layer_colonias = "colonias_cdmx"
if not proyecto.mapLayersByName(nombre_layer_colonias):
    layer_colonias = QgsVectorLayer(f"{dir_data}/colonias_iecm_2019/mgpc_2019.shp", nombre_layer_colonias)
    proyecto.addMapLayer(layer_colonias)

params_validacion = {
    'ERROR_OUTPUT' : 'TEMPORARY_OUTPUT',
    'IGNORE_RING_SELF_INTERSECTION' : False,
    'INPUT_LAYER' : layer_colonias,
    'INVALID_OUTPUT' : 'TEMPORARY_OUTPUT',
    'METHD' : 2,
    'VALID_OUTPUT' :'TEMPORARY_OUTPUT' 
}

out_validacion = processing.run("qgis:checkvalidity", params_validacion)

if out_validacion['INVALID_COUNT']:
    layer_error_existente = proyecto.mapLayersByName(out_validacion['ERROR_OUTPUT'].name())
    if layer_error_existente:
        proyecto.removeMapLayer(layer_error_existente[0])
    print('Puntos erróneos', out_validacion['INVALID_COUNT'])
    proyecto.addMapLayer(out_validacion['ERROR_OUTPUT'])

El resultado se ve así:

Se añadió al proyecto la capa (temporal) "Salida Errónea", que tiene 39 puntos señalando dónde hay errores en la capa de "colonias_cdmx". Como son muchos, no los vamos a corregir ahora.

Pero aún queremos ver el resultado de la operación de centroides a pesar de los errores, así que vamos a usar una opción que nos permite saltarse la validación y ejecutar la operación. Para eso vamos aimportar dos funciones más: QgsProcessingFeatureSourceDefinition y QgsFeatureRequest.

Creamos un nuevo script con el nombre intro_pyqgis_04.py y ponemos las siguientes líneas:

# intro_pyqgis_04.py
from qgis.core import QgsProject, QgsVectorLayer, QgsProcessingFeatureSourceDefinition, QgsFeatureRequest, QgsVectorFileWriter
import processing

proyecto = QgsProject.instance()

dir_data = "datos"
nombre_layer_colonias = "colonias_cdmx"
if not proyecto.mapLayersByName(nombre_layer_colonias):
    layer_colonias = QgsVectorLayer(f"{dir_data}/colonias_iecm_2019/mgpc_2019.shp", nombre_layer_colonias)
    proyecto.addMapLayer(layer_colonias)

params = {
    'INPUT' : QgsProcessingFeatureSourceDefinition(nombre_layer_colonias, selectedFeaturesOnly=False, featureLimit=-1, flags=QgsProcessingFeatureSourceDefinition.FlagOverrideDefaultGeometryCheck, geometryCheck=QgsFeatureRequest.GeometrySkipInvalid),
    'ALL_PARTS' : False,
    'OUTPUT' :'TEMPORARY_OUTPUT' 
}

out = processing.run("native:centroids", params)
layer_centroides = out['OUTPUT']
if not proyecto.mapLayersByName(layer_centroides.name()):
    proyecto.addMapLayer(layer_centroides)

¡Ahora sí se completó el proceso! y así se ve el resultado:

Por último, si quiero guardar el resultado de la capa de centroides puedo exportarla a un archivo con la siguiente línea:

QgsVectorFileWriter.writeAsVectorFormat(layer_centroides, f"{dir_data}/centroides_colonias.geojson",  'utf-8', layer_centroides.crs(), 'GeoJson')

Se puede exportar a varios formatos populares como shapefile, GeoJSON o GeoPackage.

Unión espacial

Por último voy a dejar otro ejemplo de un procesamiento espacial, en este caso, una unión espacial (spatial join) entre la capa de colonias y una con la ubicación de los mercados públicos de la CDMX. El resultado de esta unión es una capa temporal (Capa Unida) que contiene solo las colonias donde hay al menos un mercado público.

Creamos un nuevo script con el nombre intro_pyqgis_05.py y ponemos las siguientes líneas:

# intro_pyqgis_05.py
from qgis.core import QgsProject, QgsVectorLayer, QgsProcessingFeatureSourceDefinition, QgsFeatureRequest, QgsVectorFileWriter
import processing

proyecto = QgsProject.instance()

dir_data = "datos"
nombre_layer_colonias = "colonias_cdmx"
if not proyecto.mapLayersByName(nombre_layer_colonias):
    layer_colonias = QgsVectorLayer(f"{dir_data}/colonias_iecm_2019/mgpc_2019.shp", nombre_layer_colonias)
    proyecto.addMapLayer(layer_colonias)

nombre_layer_mercados = "mercados"
if not proyecto.mapLayersByName(nombre_layer_mercados):
    layer_mercados = QgsVectorLayer(f"{dir_data}/mercados/mercados_publicos.geojson", nombre_layer_mercados)
    proyecto.addMapLayer(layer_mercados)

params_sjoin = {
    'INPUT': QgsProcessingFeatureSourceDefinition(nombre_layer_colonias, selectedFeaturesOnly=False, featureLimit=-1, flags=QgsProcessingFeatureSourceDefinition.FlagOverrideDefaultGeometryCheck, geometryCheck=QgsFeatureRequest.GeometrySkipInvalid),
    'JOIN': layer_mercados,
    'PREDICATE': 1,
    'OUTPUT': 'TEMPORARY_OUTPUT', 
    'DISCARD_NONMATCHING': True,
}
out_sjoin = processing.run("qgis:joinattributesbylocation", params_sjoin)
layer_union = out_sjoin['OUTPUT']

layer_sjoin_existente = proyecto.mapLayersByName(layer_union.name())
if layer_sjoin_existente:
    proyecto.removeMapLayer(layer_sjoin_existente[0])

print('Features resultantes', out_sjoin['JOINED_COUNT'])
proyecto.addMapLayer(layer_union)

El resultado es el siguiente:

Conclusión

Como hemos visto, PyQGIS es una herramienta muy poderosa que nos puede ayudar a automatizar procesos al trabajar con QGIS. Puede que la curva de aprendizaje no sea la más fácil, pero con la documentación y preguntas en internet se puede llegar a realizar la mayoría de tareas cotidianas de un analista de datos geográficos.

Por último les dejo algunos recursos que he encontrado y me han parecido muy valiosos al seguir este camino:

Revisa otras entradas de este blog: