Automatizar geoprocesos con Python y GDAL/OGR
GDAL/OGR es, resumidamente, una librería de software de código abierto para operar con datos espaciales. Está presente, entre otros lugares, en sistemas de información geográfica como QGIS y permite trabajar tanto con datos ráster (GDAL) como vectoriales (OGR).
Entre las muchas formas de trabajar con GDAL/OGR se encuentra el Shell, una consola de comandos incluida en la instalación del paquete de OSGEO con la que podemos hacer toda clase de operaciones con datos espaciales a través de instrucciones. Con ella es posible trabajar con grandes volúmenes de datos que de hacerlo mediante programas con interfaz gráfica resultaría inviable.
Pues bien, una de las posibilidades que ofrece es la de ejecutar scripts de Python que automaticen la consecución de comandos necesaria para lograr nuestro objetivo. Para ello, existe una función llamada os.system() en la que podremos escribir comandos en formato texto o array para que sean ejecutados por el Shell.
En este post voy a explicar el proceso para programar un script de Python que:
El usuario tendrá que especificar al ejecutar el script:
⚠ Para que funcione correctamente, hay que ejecutar este script desde el Shell de OSGeo con python3:
C:\> python3 script.py
Si python3 no estuviera configurado y saltara error habría que escribir py3_env y ejecutar. ¡Ya estaría listo!
Al principio del script pondremos todas aquellas cuestiones necesarias para que éste funcione correctamente, como la codificación, los módulos o las variables introducidas por el usuario, así como la carpeta de salida en la que se guardará el geopackage:
# -*- coding: utf8 -*-
import os, shutil
# Valores de entrada del usuario
img_dir = str(input('Introduzca la ruta en la que se encuentran las imagenes: '))
prov_dir = str(input('Introduzca la ruta en la que se encuentran las provincias: '))
codigo = str(input('Introduzca el codigo EPSG al que quiere reproyectar las capas: '))
distancia = str(input('Introduzca la distancia entre curvas de nivel: '))
# Crear carpeta de salida
carpeta = f'{prov_dir}/cn{distancia}'
print(f'Creando la carpteta de salida -{carpeta}...')
if os.path.exists(carpeta):
shutil.rmtree(carpeta)
os.mkdir(carpeta)
print('La carpeta ya esxistia y se ha sistituido por una nueva y vacia')
else:
os.mkdir(carpeta)
print('Carpeta creada')
El siguiente paso es generar un mosaico que una todos estos MDE en uno solo, del cual extraeremos las curvas de nivel. Para ello, crearemos un ráster virtual que nos permita trabajar más adelante con todas estas imágenes pero sin generar un nuevo y pesado archivo.
La función de GDAL que lo permite es gdalbuildvrt:
print('Generando el raster virtual...')
# cambiar al directorio de las imágenes
os.chdir(img_dir)
# comando para crear un archivo de texto que liste los MDE
os.system('dir /b *.tif > archivos.txt')
# comando para generar el mosaico en forma de raster virtual
os.system('gdalbuildvrt mosaico.vrt -input_file_list archivos.txt')
print('Raster virtual generado')
A partir del ráster virtual creado anteriormente podemos extraer las curvas de nivel con gdal_contour, además de reproyectarlas con ogr2ogr al SRC que introdujo el usuario:
# Vectorizacion de curvas de nivel
print('Vectorizando curvas de nivel...')
os.system(f'gdal_contour -a elev -i {distancia} mosaico.vrt cn{distancia}.shp')
print('Curvas de nivel vectorizadas')
# Reproyeccion
print(f'Reproyectando las curvas de nivel a EPSG:{codigo}')
os.system(f'ogr2ogr -t_srs EPSG:{codigo} cn{distancia}_rep.shp cn{distancia}.shp')
print(f'Las curvas se han reproyectado a EPSG:{codigo}')
Antes de llevar a cabo la última tarea, conviene crear un índice espacial para la capa de curvas de nivel, ya que mejorará mucho el rendimiento del recorte:
print('Creando indice espacial...')
os.system(f'ogrinfo cn{distancia}_rep.shp -sql "CREATE SPATIAL INDEX ON cn{distancia}_rep"')
print('Indice espacial creado')
Ahora sí, con el siguiente código lograremos separar las curvas de nivel en capas separadas según la provincia a la que pertenezcan.
El recorte se hará también con el comando ogr2ogr
El GeoPackage es creado en el momento en el que se lleva a cabo el primer recorte. Para evitar que sea sobrescrito por los siguientes recortes y se vayan almacenando uno a uno, tendremos que hacer que se le añada automáticamente la cláusula -update una vez se haya hecho ese primer recorte. Esto lo haremos con los contadores:
print('Recortando las curvas de nivel para cada provincia en un nuevo archivo GeoPackage...')
# Cambiar a la carpeta de provincias
os.chdir(prov_dir)
# Creacion de contadores
count = 0
update = ''
# Iterar sobre los archivos de la carpeta de provincias
for provincia in os.listdir(prov_dir):
# si los archivos son shapes...
if provincia.endswith('.shp'):
# añadir -update al comando
if count > 1:
update = '-update'
nombre = os.path.splitext(provincia)[0]
print(f'Recortando curvas de nivel para la provincia de {nombre}')
# comando para recortar capas
os.system(f'ogr2ogr -f GPKG {carpeta}/cn{distancia}.gpkg -progress -clipsrc {provincia} {img_dir}/cn{distancia}_rep.shp {update} -nln {nombre}')
# incrementar el contador
count += 1
print(f'Curvas de nivel para la provincia de {nombre} recortadas')
print('Terminado')
El resultado es satisfactorio salvo por errores de geometría en la provincia de Álava y en la isla de Mallorca, pero nada grave. El GeoPackage contiene las curvas de nivel separadas según la provincia:
Tenéis el código completo en mi GitHub