-
Notifications
You must be signed in to change notification settings - Fork 0
Automatizando la exportación de mapas mediante PyQGIS
Es muy posible que, en ocasiones, os enfrentéis a la tarea siguiente: Producir mapas en los que unas capas se mantienen fijas en todos los mapas producidos y una o más capas van cambiando en cada uno de los mapas.
En mi caso tengo que producir mapas de diferentes parámetros relacionados con la investigación pesquera por especie: distribución de NASC por especie, talla media por especie, polígonos de distribución de efectivos por especie, etc. De ahora en adelante, me referiré en muchas ocasiones a la especie pero en tu caso puede ser otra magnitud o característica de los datos, como por ejemplo manzana o provincia.
Para producir cada mapa para cada especie, tenemos que activar/desactivar la capa correspondiente a la especie de interés para que se vea en la composición. El problema es que el cambio de capa no implica que la leyenda del mapa muestre la información apropiada para cada mapa. En consecuencia, es preciso retocar la leyenda cada vez que cambiamos de especie a mapear.
En ocasiones puede ser tan simple como cambiar el nombre de la capa en la leyenda. Otras veces es más complicado, ya que por cambios en los datos alteran la leyenda moviéndola, cambiando su tamaño, etc. Esto ocurre, por ejemplo, cuando usamos una simbología categorizada o graduada en la cual algún nivel desaparece.
Se puede conseguir cierta automatización o ahorro de clicks creando una composición para cada especie, por ejemplo. Pero es inevitable tener que hacer algún ajuste que otro, o tener que crear nuevas composiciones si aparecen especies nuevas, y esto suele llevar a errores (Yo, sin ir mas lejos, produzco unos 70 mapas de este tipo al año y todos los años tengo que revisar alguno por pequeños defectos)... y aun con cierto ahorro de tiempo, ¡tienes que generar todos los mapas a mano!
Aquí es donde PyQGIS nos echa un cable: Podemos crear scripts en Python que automaticen toda esta tarea y librarnos de mucho trabajo repetitivo.
- Las instrucciones que voy a dar han sido testadas con QGIS 3.4.
- PyQGIS nos permite realizar estas tareas tanto mediante la consola Python de QGIS como en scripts/aplicaciones standalone (independientes), lo que hace posible crear aplicaciones QGIS fuera de QGIS, por así decirlo. Las instrucciones que os daré a continuación están pensadas para hacerlo en modo independiente, pero son fácilmente adaptables a la consola QGIS. De hecho, la única diferencia se encuentra en la importación de las librerías al comienzo y en la forma de acceder al proyecto.
- Mediante PyQGIS podríamos crear nuestro proyecto y composiciones de mapa desce cero, ajustando todos los parámetros que quisiéramos, cargando y eliminando capas, etc., programáticamente. Pero nosotros vamos a aprovechar las ventajas que nos ofrece contar con un GUI y crearemos una proyecto y unas composiciones de mapa que nos servirá de plantilla. De esa manera nos ahorraremos un montón de ajustes y de código.
Bueno, pues lo primero que tenemos que hacer es cargar el módulo de QGIS e iniciar la aplicación QGIS. Hay varias maneras de hacerlo, que podéis consultar aquí
Pero quizá la más popular, y es la que voy a emplear, es esta:
from qgis.core import *
QgsApplication.setPrefixPath("", True)
qgs=QgsApplication([],False)
qgs.initQGis()
Este código no es necesario si lo haces dentro de la consola Python de QGIS, ya la consola Python importa estas librerías al lanzarse.
Ahora tenemos que acceder al proyecto. Esta es la otra parte que es diferente en la aplicación independiente y desde la consola. Para acceder al proyecto desde la consola de QGIS es muy sencillo, simplemente reclamamos la instancia del proyecto:
proyecto=QgsProject.instance()
En el caso de nuestra aplicación independiente, es un poquito más complicado. Tenemos que leer el proyecto desde un archivo, pero primero tenemos que crear el objeto proyecto y después leer el proyecto:
proyecto=QgsProject()
proyecto.read('/home/fulano/GIS/miproyecto.qgz')
A partir de ahora, el código es idéntico tanto para la consola de QGIS como para una hipotética aplicación independiente.
El siguiente paso es acceder al gestor de composiciones de mapa de nuestro proyecto.
composiciones=proyecto.layoutManager()
Continuamos accediendo a la composición que nos interese, que previamente habremos diseñado mediante el compositor de mapas de QGIS y guardado en el proyecto que hemos cargado.
En este ejemplo, aparte de una capa con las masas de tierra y otra de batimetría, disponemos de una leyenda y hemos proporcionado un marco, una escala y, en definitiva, hemos creado una composición a nuestro gusto.
composicion=composiciones.layoutByName('PRUEBA')
Bien... lo que más nos interesa modificar para cada mapa es la leyenda. QGIS guarda los elementos de la composición (leyendas, rejillas, mapas, escalas) en una lista. Desde la consola es fácil distinguir la leyenda de los demás ya que podemos hacer:
items=composicion.items()
examinar la clase de los items y quedarnos con el que corresponda a un objeto de la clase QgsLayoutItemLegend, si resultara ser el elemento 3 de la lista:
leyenda=composicion.items[3]
No obstante, si hacemos una aplicación independiente no es tan sencillo ya que no podemos examinar la lista. En este caso, dado que lo más habitual es que la composición tenga una única leyenda, haremos:
leyenda=[item for item in items if type(item)==QgsLayoutItemLegend][0]
Para cualquier otro elempento de la composición podermos hacer lo mismo. Por ejemplo, si nos interesa acceder al mapa, en el caso de que solo sea uno, podemos hacer análogamente:
mapa=[item for item in items if type(item)==QgsLayoutItemMap][0]
Aprovecho para decir que acceder al mapa es fundamental para poder colocar bien la leyenda posteriormente .
Como hemos dicho antes, lo suyo es crear el proyecto y la composición en QGIS y luego aprovecharlo en la apicación. Al crear la composición de mapa nos conviene que la leyenda cuente con todo el formateo posible (tipos de letra del título, centrado del mismo, borde, transparencias) ya que esto es mucho más complicado cambiarlo programáticamente.
QGIS almacena la leyenda y las capas de la vista de mapa como árboles mediante un objeto de tipo QgsLayoutItemLegend que tiene un modelo asociado.
Hay un efecto secundario de trabajar sobre la leyenda y es que ciertos cambios en la leyenda se propagan hacia el arbol de capas del proyecto. Por eso es muy importante impedir que la leyenda se autoactualice, porque en ese caso al hacer la operación clear() en el árbol de capas de la leyenda del mapa nos cargamos el árbol de capas del proyecto.
Por lo tanto, lo primero que haremos antes de modificar nada de la leyenda es:
leyenda.setAutoUpdateModel(False)
Ahora accedemos al modelo de la leyenda y, de ahí, al árbol en sí:
modelo=leyenda.model()
arbol=modelo.rootGroup()
Es conveniente limpiar el árbol de la leyenda. Recuerda que el resto de elementos y propiedades (título,marco, transparencia, tipos de letra...) se mantienen
arbol.clear()
Bien, en mi caso la leyenda consta de título (que ya lo hemos generado) y de un grupo con la capa de información por especie y otro con la capa de batimetría. Crearemos esos grupos asignándolos a una variable para trabajar mejor después:
nascgrupo=arbol.addGroup('NASC, m²/nm²')
batigrupo=arbol.addGroup('Bathymetry')
Ahora localizamos las capas que sea necesario para añadirlas a cada grupo. Todo esto es mucho más fácil si se garantizan nombres únicos de capas. Para localizar una capa, si es que tiene el nombre único:
proyecto.mapLayersByName('NASC_ANE_ECO2017')
Esto nos devuelve una lista con los nombres de capa que coincidan con el criterio de búsqueda. Como se supone, en el caso de que las capas tengan nombre único, que sólo debe haber una, podemos hacer
capa=proyecto.mapLayersByName('NASC_ANE_ECO2017')[0]
Si hemos sido espabilados y tenemos todos los nombres de capa en una lista, podemos obtener rápidamente las capas:
nombresCapa=['NASC_ANE_ECO2017', 'NASC_PIL_ECO2017', 'NASC_HOM_ECO2017']
capasMapa=[proyecto.mapLayersByName(nombre)[0] for nombre in nombresCapa]
Esta lista de capas será sobre la que iteremos posteriormente para obtener nuestros mapas por especie (o, te recuerdo, por manzana, provincia o lo que sea)
También tenemos que recuperar aquellas capas que queremos que aparezcan en la leyenda pero que no sean sobre las que iteramos, como por ejemplo batimetría, que es fija en todos los mapas y que ha de aparecer en la leyenda para que se vea el color correspondiente a cada profundidad.
baticapa=proyecto.mapLayersByName('BATIMETRIA_GENERAL_GC')[0]
Para hacer visible o no una capa, es un poco complicado pero muy parecido a la manera de trabajar con la leyenda. Se tiene que acceder al arbol de capas del proyecto, y localizar la capa, lo que ocurre es que en este caso es más conveniente localizar la capa por ID. Para eso nos valemos del método id() de la capa o, directamente, pasando la capa a findLayer()
proyecto.layerTreeRoot().findLayer(<capa>).setItemVisibilityChecked(True)
Lo suyo es en principio dejar visibles sólo las capas de cobertura de interés y luego ir iterando con el resto.
Para incorporar o quitar una capa de la leyenda del mapa, sólo tenemos que usar el método apropiado del grupo de capas al quenos interesa añadir la capa. Recordemos que anteriormente creamos el grupo batigrupo así que, por ejemplo, añadimos una capa al grupo mediante:
batigrupo.addLayer(baticapa)
Y la eliminamos con:
batigrupo.removeLayer(baticapa)
Cuando hacemos addLayer estamos añadiendo un nodo al arbol de capas y grupos de la leyenda. Por eso es preferible asignar la salida de addLayer a una variable, así es mucho más fácil cambiar cosas en caso de necesidad, como por ejemplo:
nodobati=batigrupo.addLayer(baticapa)
nodobati.setName('')
Haciendo el nombre de la capa='' se consigue que no se vea en la leyenda. El problema es que se modifica también en la lista de capas del proyecto (no sé muy bien por qué). Esto se puede conseguir también de esta manera, no tan intuitiva:
QgsLegendRenderer.setNodeLegendStyle(nodobati,QgsLegendStyle.Hidden)
Ahora, con nuestras capas y grupos en su sitio, tenemos que generar los mapas. El proceso, en general, es sencillo e intuitivo:
- Hacemos visible una capa
- Generamos el mapa correspondiente
- Hacemos invisible esa capa
- Hacemos visible la siguiente capa
- ....
Empezamos creando una lista con las capas. Si tenemos los nombres de las capas y éstos son únicos, podríamos hacer como describimos anteriormente:
nombresCapa=['NASC_ANE_ECO201907', 'NASC_PIL_ECO201907',
'NASC_MAS_ECO201907', 'NASC_MAC_ECO201907', 'NASC_HOM_ECO201907',
'NASC_HMM_ECO201907']
capasMapa=[]
for nombre in nombresCapa:
capasMapa.append(proy.mapLayersByName(nombre)[0]
o también:
capasMapa=[proy.mapLayersByName(nombre)[0] for nombre in nombresCapa]
Para producir la salida, lo suyo es , en primer lugar, quitar la visibilidad de todas las capas de interés
for capa in capasMapa:
proy.layerTreeRoot().findLayer(capa).setItemVisibilityChecked(False)
Ahora tenemos que crear el exportador de la composición de mapas, que es quien se encargará de hacer la magia que genera los pdf, SVG o imágenes:
exporter=QgsLayoutExporter(composicion)
Uno de los mayores problemas que provoca tener que hacer muchos ajustes manuales cuando exportamos los mapas desde QGIS es, como dije anteriormente, el cambio de tamaño de la leyenda, que hace que se desajuste de su posición. Para corregirlo en nuestro script, necesitamos saber la posicion del mapa en la página para poder establecer bien la posicion de la leyenda. Queremos dejarla a la misma distancia del borde inferior del mapa que del borde izaquierdo Para evitar estar haciendo cálculos todo el rato, vamos a fijar el punto de referencia para los movimientos de la leyenda a la esquina inferior izquierda de la misma. De esta manera, conociendo la distancia horizontal entre mapa y leyenda daremos la vertical, pero primero hemos de asegurarnos de que el punto es el superior izquierdo para hacer los cálculos. Posteriormente nos aseguraremos de que sea el inferior izquierdo.
leyenda.setReferencePoint(QgsLayoutItem.UpperLeft)
posMapa=mapa.pagePositionWithUnits()
tamMapa=mapa.sizeWithUnits()
mapaX=posMapa.x()
mapaY=posMapa.y()
Como queremos calcular la esquina inferior izquierda, sólo necesitamos calcular la coordenada y
mapaY2=mapaY+tamMapa.height()
A partir de la posición de la leyenda sacamos el margen entre leyenda y mapa, que será la distancia a aplicar
leyendaX=leyenda.pagePositionWithUnits().x()
distanciaMargen=abs(leyendaX-mapaX)
Con todo lo anterior, sacamos el punto en el que anclar la lenda
anclaX=mapaX+distanciaMargen
anclaY=mapaY2-distanciaMargen
anclaLeyenda=QgsLayoutPoint(anclaX,anclaY)
Volvemos a hacer el punto de referncia abajo a la izquierda
leyenda.setReferencePoint(QgsLayoutItem.LowerLeft)
Y ya con esto, iteramos con las capas
for capaTrabajo in capasMapa:
#La hacemos visible en el proyecto y todos sus ancestros, para que se vea
proyecto.layerTreeRoot().findLayer(capaTrabajo).setItemVisibilityCheckedParentRecursive(True)
nascgrupo.addLayer(capaTrabajo) # La metemos en la leyenda
# Creamos el nombre de fichero a partir del nombre de la capa
nombreImagen='/home/fulano/GIS/MAPA_%s.png' %capaTrabajo.name()
leyenda.adjustBoxSize()
leyenda.refresh()
# Ahora reposicionamos la leyenda
leyenda.attemptMove(anclaLeyenda, useReferencePoint=True, includesFrame=True)
leyenda.positionAtReferencePoint(6)
# Y producimos la salida a imagen
exporter.exportToImage(nombreImagen, exporter.ImageExportSettings())
# Ahora hacemos la capa de trabajo invisible, pero ojo, no hacemos invisibles sus ancestros ya que
# eso incluye a las capas fijas (batimetria, etc)
proy.layerTreeRoot().findLayer(capaTrabajo).setItemVisibilityChecked(False)
# Y finalmente sacamos la capa de la leyenda
nascgrupo.removeLayer(capaTrabajo)
A la hora de iterar existen otras aproximaciones, como podría ser acceder a un grupo concreto de capas del proyecto e iterar directamente sobre las capas de ese grupo.
Sí que os diré que la parte del posicionamiento funciona de una manera un poco errática y me ha costado bastante dar con la manera correcta de hacerlo e ignoro si es así de complicado o se debe a algún problema adicional.
Y ahora viene el regalito...
Bien, esta aplicación podría considerarse, sin lugar a dudas, una espantajería o una gran idea. Yo no soy programador y quizá lo que os propongo ahora es aberrante o hay mejor manera de hacerlo, seguro. Pero mi interés en ella se centra en que aunque R produce unas salidas en ocasiones espectaculares en materia de mapas, crear mapas no es tan sencillo y requiere de bastante dominio de R. Crear mapas como los que se pueden crear con QGIS con sólo unos clicks y añadiendo rejillas diversas, escalas, leyendas variopintas... es algo que no digo que no se pueda con R pero desde luego no está al alcance de todos...
Pero, por otro lado, la existencia de RMarkdown hace muy sencillo y versátil la producción de informes. Además, se puede mezclar con LaTex y en definitiva es una herramienta tremenda para esto.
Incorporar nuestros mapas de QGIS a informes RMarkdown es muy sencillo. Aquí no me refiero a cargar las imágenes desde archivos que hayamos generado antes, sino a incorporar en el documento RMarkdown nuestro script de PyQgis y generar los mapas al vuelo.
Para ello haremos uso en R de la librería Reticulate. Este ejemplo que os pongo a continuación es experimental y tengo que segir investigando, pero es muy prometedor y funciona ¡Probad a ver!
---
title: "Untitled"
output: pdf_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(reticulate)
use_python('/usr/bin/python3')
```
```{python, engine.path = '/usr/bin/python3'}
from qgis.core import *
QgsApplication.setPrefixPath("", True)
qgs=QgsApplication([], False)
qgs.initQgis()
pro=QgsProject()
pro.read('myproject.qgs')
la=pro.layoutManager().layoutByName('MYLAYOUT')
exporter=QgsLayoutExporter(la)
img=exporter.renderPageToImage(0)
img.save('./mymap.png')
```
\includegraphics{mymap.png}
Bueno, pues hasta aquí hemos llegado. Un saludo para todos de @imasdemase