Por Rober J
Python y ArcGIS
Cada SIG cuenta con su propia librería de Python que permite acceder a los geoprocesos de dicho SIG. En este caso, ArcPy es la librería que da acceso a las funciones de ArcGIS en un entorno Python, dándonos acceso a las cajas de herramientas de geoprocesamiento estándar y a la posibilidad de usar otros módulos (siempre que tengamos licencia para usarlos)
¿Qué es lo que cambia? Parece contraintuitivo sustituir un amigable cuadro de texto por un churro de texto, pero gracias a ello accedemos a una herramienta mucho más flexible, ya que dentro de un script (un pequeño código) podemos diseñar qué se ejecuta, cuándo se ejecuta y con qué parámetros, encadenando unos procesos con otros y obteniendo resultados a nuestra medida.
Por ejemplo, ese ‘output’ que hemos especificado en la función arcpy.Clip podemos meterlo a continuación en otra función distinta, o en varias, y esos resultados pasarlos por otro geoproceso y así sucesivamente.
Cómo comenzar a usar ArcPy
Puede hacerse de varias maneras, principalmente:
A través de la ventana de Python incorporada en ArcMap. Es una opción rápida para ejecutar pequeños scripts sin complicarnos demasiado, pero se queda corto puesto que su funcionalidad se reduce a la de escribir código o pegarlo desde otra fuente y ejecutarlo.
En este caso he ejecutado el Clip de la imagen anterior, recortando la red fluvial usando el polígono de un municipio cualquiera:
Haciendo uso de un entorno de desarrollo integrado (IDE), es decir, un software diseñado para trabajar con código y que dispone de herramientas para hacernos la vida más fácil: predicción de texto, resalte de sintaxis, documentación instantánea, cambios múltiples…
Existen múltiples IDEs y es cosa de cada uno escoger el que le vaya bien. Personalmente he usado Microsoft Visual Studio Code y PyCharm, y me quedo con el segundo puesto que está diseñado específicamente para Python y me ha dado menos problemas a nivel general (todos tienen sus cosillas…)
Aspecto de PyCharm
En cualquier caso, deberás configurar el IDE para que tenga acceso a los módulos de ArcGIS vinculándolo con el intérprete que ArcGIS trae por defecto. El intérprete es el programa que traduce el código de Python para que el ordenador pueda ejecutarlo, y el módulo ArcPy solo funcionará si se utiliza el IDE junto a este intérprete.
Generalmente, ArcGIS instala Python en una carpeta llamada Python27 en la raíz del disco duro (generalmente suele ser C:\\ ) por lo que habrá que buscar en esa carpeta el archivo python.exe y seleccionarlo como intérprete.
Obtener scripts de Model Builder
Una forma de comenzar montar un código de Python para ArcGIS es utilizar Model Builder para conceptualizar el trabajo que queremos hacer y extraer de él las funciones de geoproceso que necesitaremos.
En el artículo Análisis de ubicación de un vertedero con Model Builder comento brevemente las ventajas y limitaciones que presenta Model Builder y sus semejanzas con Python. Pues bien, podemos exportar los modelos a archivos Python (archivos con extensión .py) y abrirlos con una IDE para editarlos.
Siguiendo con el ejemplo del Clip, he construido el modelo en Model Builder y lo he exportado de la siguiente manera:
A continuación lo abro en un IDE y presenta el siguiente aspecto. Como comentaba, es muy interesante porque te da una estructura básica a partir de la cual continuar desarrollando el script, pero habrá que modificarlo para que funcione:
# -*- coding: utf-8 -*- # --------------------------------------------------------------------------- # clip.py # Created on: 2021-03-04 22:32:42.00000 # (generated by ArcGIS/ModelBuilder) # Description: # --------------------------------------------------------------------------- # Import arcpy module import arcpy # Local variables: Red_fluvial = "Red_fluvial" seleccion = "seleccion" output = "C:\\Users\\Roberto\\Documents\\ArcGIS\\Default.gdb\\output" # Process: Clip arcpy.Clip_analysis(Red_fluvial, seleccion, output, "")
Estructura de un script
Del código anterior podemos diferenciar varias partes que funcionan a modo de esqueleto para un script:
- La codificación de caracteres que va a usarse (utf-8)
- Metadatos del archivo .py como su nombre o fecha de creación
- Importación del módulo arcpy (¡no olvidar!)
- Variables locales: conjunto de variables que definen los parámetros de los geoprocesos, como son las rutas de las capas de entrada (Red_fluvial y seleccion) y de salida (output)
- Procesos: la parte del código que ejecutará las funciones de geoproceso haciendo uso de las variables que definimos en el punto anterior
Sin embargo, esta estructura base no funcionará correctamente fuera del entorno de ArcMap porque no reconocerá las capas, por lo que tendremos que modificar las variables añadiendo una ruta válida. Además, aunque cambiemos la ruta, este script solo podrá ejecutarse 1 sola vez porque se generaría una capa con un nombre que ya existe, por lo que tendremos que asegurarnos de tener correctamente configurados algunos parámetros de variables de entorno como veremos a continuación.
Importación de módulos
Antes de nada, tendremos que importar ArcPY junto al resto de módulos que vayamos a usar durante el script. Si, por ejemplo, vamos a querer que se cree automáticamente una nueva carpeta con el resultado, no se nos puede olvidar importar también el módulo os:
import arcpy, os
Importar arcpy tal como mostramos arriba da acceso a las siguientes funcionalidades (fuente: Esri):
- Caja de herramientas de Análisis (Analysis Tools)
- Caja de herramientas Cartografía (Cartography tools)
- Caja de herramientas Conversión (Conversion tools)
- Caja de herramientas Administración de datos (Data Management Tools)
- Caja de herramientas Edición (Editing Tools)
- Caja de herramientas Geocodificación (Geocoding Tools)
- Caja de herramientas Referencia lineal (Linear Referencing Tools)
- Caja de herramientas de Multidimensión (Multidimension Tools)
- Caja de herramientas Estadística espacial (Spatial Statistics Analyst)
Sin embargo, existen otros módulos que deben importarse a parte para acceder a más funciones de ArcGIS:
- arcpy.da – módulo de acceso de datos
- arcpy.mapping – módulo de representación cartográfica
- arcpy.sa – módulo del Spatial Analyst
- arcpy.na – módulo del Network Analyst
- arcpy.ga – módulo del Geostatistical Analyst
Para llevar a cabo las operaciones con geodatos vectoriales que presento a continuación se deben importar los siguientes módulos y definir las siguientes variables de entorno:
# Modulos import arcpy from arcpy import env # Entorno ruta = 'C:\\...' env.workspace = ruta env.overwriteOutput = True
Para llevar a cabo las operaciones con datos ráster que presento a continuación se deben importar los siguientes módulos y definir las siguientes variables de entorno:
# Modulos
import arcpy
from arcpy import env
from arcpy.sa import *
# Entorno
ruta = 'C:\\...'
env.workspace = ruta
env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")
Variables de entorno
Las variables de entorno o environments son unos parámetros o funciones que conviene definir al comienzo del script (justo tras la importación de módulos) para que los geoprocesos funcionen de una u otra manera. Son, por así decirlo, ‘las reglas’ que regirán el script.
Estas variables se encuentran dentro de la clase env de ArcPy. Son bastante numerosas y podéis encontrarlas todas aquí, pero los más básicos serían:
- Directorio de trabajo – la carpeta en la que se localizan los geodatos. Definirlo es útil porque nos permitirá ahorrarnos el tener que escribir rutas completas en el futuro, es decir, podremos llamar a las capas solo por su nombre y su extensión, ya sean inputs u outputs.
- Sobreescritura de archivos – al definirla como True se borrarán de forma automática las capas antiguas que tengan el mismo nombre que una capa nueva que se acabe de generar. En nuestro ejemplo del Clip, al tener este parámetro activado el segundo Clip borraría el primero ya que el output en este caso tiene siempre el mismo nombre.
- Sistema de proyección – establecer el SRC de nuestro marco de datos. Al igual que en ArcMap, se proyectarán las capas ‘al vuelo’ usando el SRC de la primera capa leída por nuestro script.
- Activación de extensiones – muchos geoprocesos como los del módulo Spatial Analyst se encuentran bajo licencia, por lo que deben activarse del mismo modo que hacemos en ArcMap – Customize – Extensions…
# Definir el directorio de trabajo
arcpy.env.workspace = 'ruta'
# Activar la sobreescritura de archivos
arcpy.env.overwriteOutput = True
# Establecer el SRC al ETRS89 UTM Zona 30 Norte
arcpy.env.cartographicCoordinateSystem = "Coordinate Systems\Projected Coordinate Systems\UTM\Europe\ETRS 1989 UTM Zone 30N.prj"
# Activar la extensión Spatial Analyst
arcpy.CheckOutExtension('spatial')
Entre otros entornos están el de establecer un SRC para las capas de salida, la resolución de las nuevas capas ráster, crear pirámides o el añadir las nuevas capas a la visualización.
Variables locales
Como dijimos antes, son el conjunto de variables que usarán los geoprocesos para llevarse a cabo. Suelen ser variables locales:
- Las rutas de las capas de entrada
- Las rutas de las capas de salida
- Filtros de archivos
- Cálculos de valores
- Expresiones SQL para hacer selecciones
- …
Hay tantas variables variables locales como distintos geoprocesos que vayamos a utilizar y parámetros de éstos tengamos que introducir. Para la definición de estas variables, es habitual usar inputs para que sea el usuario el que las defina sobre la marcha.
Geoprocesos
Es al final del script cuando deberíamos colocar los geoprocesos, puesto que éstos harán uso de las variables que hemos definido previamente.
La sintaxis de todos estos procesos se encuentran en la documentación oficial de Esri a la que podéis acceder a través de los enlaces del apartado Importación de módulos.
Algunas de las funciones de geoproceso más básicas de ArcPy son:
## Crear una capa temporal a partir de una capa existente
arcpy.MakeFeatureLayer_management('capa_entrada', 'lyr')
## Seleccionar entidades según sus atributos
arcpy.SelectLayerByAttribute_management('lyr', 'TIPO_SELECCION', 'expresión SQL')
## Copiar entidades a una nueva capa
arcpy.CopyFeatures_management('lyr', 'nueva_capa')
## CLIP
arcpy.Clip_analysis('capa_entrada', 'capa_clip', 'capa_salida', ' ')
## BUFFER
arcpy.Buffer_analysis('capa_entrada', 'capa_salida', 'distancia', 'FULL', 'ROUND', 'NONE', ' ', 'PLANAR')
## CALCULATE STATISTICS (raster)
arcpy.CalculateStatistics_management('capa_entrada')
## RESAMPLE (modificar resolución)
arcpy.Resample_management('capa_entrada', 'capa_salida', 'nueva_resolución', 'MÉTODO')
Comprobar si un campo existe, y si no crearlo
Para ello hay que usar la función arcpy.ListFields() para obtener una lista con objetos de tipo campo correspondientes a los campos de la tabla de atributos de una capa y comprobar su existencia. Si no existe, usar la función arcpy.AddField_management() para agregar el campo
nuevo_campo = 'nombre_campo'
capa = 'nombre_capa.shp'
listaCampos = arcpy.ListFields(capa)
existencia = 0
for campo in listaCampos:
if campo.name == nuevo_campo:
existencia = 1
if existencia == 1:
print('El campo ' + nuevo_campo + ' ya existe')
else:
arcpy.AddField_management(capa,nuevo_campo,tipo...)
print('El campo ' + nuevo_campo + ' ha sido creado')
Listar capas vectoriales de un directorio
La función ListFeatureClasses crea listas con los nombres de las capas junto a su extensión que se encuentran en el directorio de trabajo definido en las variables de entorno. Además, permite filtrarlas por nombre y tipo:
# Listar todas las capas vectoriales
arcpy.ListFeatureClasses()
# Listar solo las capas de puntos
arcpy.ListFeatureClasses(,'Point')
# Listar solo las capas de líneas
arcpy.ListFeatureClasses(,'Line')
# Listar solo las capas de polígonos
arcpy.ListFeatureClasses(,'Polygon')
# Listar solo las capas cuyo nombre empiece por 'Col'
arcpy.ListFeatureClasses('Col*')
# Listar solo las capas cuyo nombre termine por por 'egios'
arcpy.ListFeatureClasses('*egios')
# Listar solo las capas cuyo nombre coincida con 'Colegios'
arcpy.ListFeatureClasses('Colegios')
# Listar solo las capas cuyo nombre coincida con 'Colegios' y sean de tipo poligonal
arcpy.ListFeatureClasses('Colegios', 'Polygon')
Crear capas temporales
Las capas temporales o capas layer que solo existen mientras se ejecuta el script y nos permiten hacer selecciones y otras operaciones sin modificar la capa original. El primer argumento es para la capa que vamos a 'duplicar' y el segundo para darle el nombre con el que se identificará durante el script:
arcpy.MakeFeatureLayer_management("capa_entrada.shp", "capa_lyr")
No es necesario indicar la extensión de la capa temporal.
Seleccionar entidades por atributos y guardarlas en una capa nueva
El siguiente código crea una nueva selección con los árboles que miden más de 15 metros y los guarda en un archivo nuevo.
arcpy.MakeFeatureLayer_management("arboles.shp", "arboles_lyr")
arcpy.SelectLayerByAttribute_management('arboles_lyr', "NEW_SELECTION", '"ALTURA" > 15')
arcpy.CopyFeatures_management('arboles_lyr', 'arboles_15m.shp'')
En la documentación de Esri tenéis más detalles sobre la función SelectLayerByAttribute
Seleccionar entidades por localización y guardarlas en una capa nueva
También se pueden hacer selecciones espaciales basadas en las distintas relaciones espaciales entre entidades geográficas. En el siguiente ejemplo se seleccionan aquellos árboles que se encuentren dentro de una determinada parcela guardada como capa individual:
arcpy.MakeFeatureLayer_management('arboles.shp', 'arboles_lyr')
SelectLayerByLocation('arboles_lyr', 'WITHIN', 'parcela.shp', 'NEW_SELECTION')
arcpy.CopyFeatures_management('arboles_lyr', 'arboles_parcela.shp'')
En la documentación de Esri tenéis más detalles sobre la función SelectLayerByLocation
Cursores - acceder a las tablas de atributos
Los cursores crean listas con objetos de tipo fila. Cada uno de estos objetos nos permite acceder a los valores de los distintos campos que contiene una fila.
Son necesarios para acceder a los registros de las tablas de atributos de nuestras capas. Hay 3 tipos:
- arcpy.SearchCursor() es el cursor de búsqueda, válido para hacer consultas sin modificar valores
- arcpy.UpdateCursor() es el cursor de actualización necesario para hacer cambios en los registros existentes
- arcpy.InsertCursor() es el cursor para insertar nuevos registros en una tabla de atributos
A continuación tenéis ejemplos del uso de los cursores:
Para hacer esto hay que aplicar el cursor de búsqueda arcpy.SearchCursor() a la capa de la queramos leer la información y guardarlo como objeto. Después habrá que llamar a este objeto para que se imprima el campo indicado dentro de un bucle:
cursor = arcpy.SearchCursor('nombre_capa.shp')
for fila in cursor:
print(fila.nombre_campo)
Para ello hay que guardar como objeto el cursor de actualización arcpy.UpdateCursor() aplicado a la capa que queremos modificar. Después deberemos pasarle a este objeto un bucle for o while.
En este caso se usa un bucle for para cambiar todos los valores de un campo por el número 3:
cursor = arcpy.UpdateCursor('nombre_capa.shp')
for fila in cursor:
fila.campo = 3
cursor.updateRow(fila)
También se podría usar en su lugar un bucle while. Para que este bucle itere sobre cada fila hay que indicarle que vaya moviéndose a la siguiente posición del cursor a cada paso que da con el método .next. Cuando llega al final no devuelve ningún valor y el bucle se cierra:
cursor = arcpy.UpdateCursor('nombre_capa.shp')
fila = cursor.next()
while fila:
fila.campo = 3
cursor.updateRow(fila)
fila = cursor.next()
Al acabar de hacer cambios se deberían eliminar los objetos creados para evitar errores en el almacenamiento de los nuevos valores:
del cursor del fila
Se puede, por ejemplo, cambiar solo algunos valores, como por ejemplo los que sean mayores que 7 para que pasen a ser 0:
cursor = arcpy.UpdateCursor('nombre_capa.shp')
for fila in cursor:
if fila.prueba > 7:
fila.prueba = 0
cursor.updateRow(fila)
else:
continue
del cursor
Además de modificar quirúrgicamente los valores que nos interesan de un campo concreto, podemos usar este procedimiento para reclasificar todos los valores de un campo, algo muy útil cuando queremos transformar variables continuas en discretas:
cursor = arcpy.UpdateCursor('estaciones_meteo.shp')
for fila in cursor:
if fila.temperatura < 10:
fila.descripcion = 'frío'
cursor.updateRow(fila)
elif fila.temperatura >= 10 AND fila.temperatura < 25:
fila.descripcion = 'templado'
cursor.updateRow(fila)
else:
fila.descripcion = 'cálido'
cursor.updaterow(fila)
del cursor
Listar y filtrar capas ráster
La función ListRasters crea listas con los nombres de las capas junto a su extensión que se encuentran en el directorio de trabajo definido en las variables de entorno. Además, permite filtrarlas por nombre y tipo:
# Listar todos los rasters
arcpy.ListRasters()
# Listar solo los rasters de tipo TIFF
arcpy.ListRasters(,'tif')
# Listar solo los rasters GRID
arcpy.ListRasters(,'grid')
# Listar solo los rasters IMG
arcpy.ListRasters(,'img')
# Listar solo los rasters cuyo nombre empiece por 'Temp'
arcpy.ListRasters('Temp*')
# Listar solo las capas cuyo nombre termine por por 'maximas'
arcpy.ListRasters('*maximas')
# Listar solo las capas cuyo nombre coincida con 'Temperaturas maximas'
arcpy.ListRasters('Temperaturas maximas')
# Listar solo las capas cuyo nombre empiece por 'Temp' y sean de tipo TIFF
arcpy.ListRasters('Temp*', 'TIFF')
Obtener la resolución de una capa
Las unidades del valor devuelto variarán en función de la proyección de la capa (grados, metros...)
capa_raster = Raster('MDT.tif')
resolucion = capa_raster.meanCellWidth
Calcular estadísticas
Este geoproceso habilita a las capas ráster para aplicar posteriormente algunas herramientas de ArcPy y evitar el error 001100: Failed because no statistics is available
arcpy.CalculateStatistics_management(capa_raster)
Obtener las coordenadas de los límites de la capa
Con la función Point se puede obtener un objeto de tipo punto con los valores devueltos por los métodos .extent:
capa_raster = Raster('MDT.tif')
limites = 'XMIN = {0}, XMAX = {1}, YMIN = {2}, YMAX = {3}'.format(capa.extent.XMin,capa.extent.XMax,capa.extent.YMin,capa.extent.YMax)
print(limites)
Este código devolverá algo así:
XMIN = 569301.0, XMAX = 810701.0, YMIN = 4413136.0, YMAX = 4755136.0
Reclasificar valores - álgebra de mapas
Para realizar operaciones sobre los pixeles de un ráster se debe usar la función Raster() sobre la imagen que se quiere reclasificar y a continuación utilizar los operadores de Python para modificar los valores.
Hay que recordar que las reclasificaciones se guardan como variables que deben ser almacenadas en archivos nuevos usando el método .save()
En el siguiente ejemplo se reclasifica un mapa de radiación solar generando una nueva imagen en la que los valores superiores a los 5 kWh/m2 pasen a tener valor 1 y el resto 0:
reclasificacion = Raster('radiacion.tif') > 5
reclasificacion.save('radiacion_reclas.tif')
También se pueden reconvertir las unidades, convirtiendo los kWh en megajulios multiplicando los valores por 3.6:
conversion = Raster('radiacion.tif') * 3.6
conversion.save('radiacion_mjulios.tif')
Además, se pueden usar varios operadores combinando consultas separadas por paréntesis. En el siguiente ejemplo se van a reclasificar los valores de un mapa de orientaciones para que los valores entre el noreste y el sureste (entre 45º y 135º) cambien a valor 1 y el resto sean 0:
orientacion_este = (Raster('aspect.tif') > 45) & (Raster('aspect.tif') < 135)
orientacion_este.save('aspect_este.tif')
Incluso se pueden combinar los datos de varias capas ráster para combinar variables distintas. Con este código pasarían a valer 1 los píxeles de un MDT que se encuentran entre los 300 y los 500 metros de altitud y que además se encuentren en una pendiente superior a los 5º tomada de otra capa ráster:
altitud_pendientes = (Raster('MDT.tif') > 300) & (Raster('MDT.tif') < 500 & (Raster('slope.tif') > 5)
altitud_pendientes.save('reclasificacion.tif')
🗺 Nota: si se usan capas distintas éstas deben encontrarse en el mismo SRC, o si no los píxeles no coincidirán
Matrices
Operar con matrices reduce el tiempo de procesado de las operaciones sobre capas ráster. Consiste en convertir la información del ráster, que es en sí una matriz de valores, en una matriz o 'tabla' con la que poder operar usando el módulo NumPy de Python (activado por defecto).
Fuente: docs.qgis
Para convertir una capa ráster en una matriz de NumPy, ArcPy posee la función arcpy.RasterToNumPyArray(). En el siguiente ejemplo se lleva a cabo dicha conversión especificando que los valores nulos del ráster tengan valor 0:
capa_raster = raster('MDT.tif')
raster_matriz = arcpy.RasterToNumPyArray(capa_raster,nodata_to_value=0)
Una vez hecho esto, podremos usar el objeto o variable raster_matriz que hemos creado para aplicarle los métodos propios del módulo NumPy. Tenéis a continuación ejemplos de operaciones usando algunos de estos métodos:
# Sumar todos los valores de lo píxeles raster_matriz.sum() # Extraer el valor mínimo raster_matriz.min() # Extraer el valor máximo raster_matriz.max()
Cambiar la resolución de una capa
El siguiente código comprueba si el tamaño del píxel de una capa ráster supera un umbral que hayamos establecido. En tal caso, se creará una nueva imagen con la resolución establecida en dicho umbral:
capa_raster = Raster('nombre_archivo.extensión')
umbral = X
resolucion = capa_raster.meanCellWidth
print('La resolución del ráster es ' + str(resolucion))
if resolucion > umbral:
print('La resolución es demasiado baja. Aumentando la resolución a' + str(umbral) + ' metros' )
arcpy.Resample_management(capa_raster, 'nombre_salida_' + str(umbral) + '.extensión', umbral, 'NEAREST')
nuevo_raster = Raster('nombre_salida_' + str(umbral) + '.extensión')
else:
print('Resolución óptima')
En la documentación de Esri tenéis más detalles sobre la función Resample_management