Operaciones sobre entidades con QGIS y Python

Cuando hay que realizar muchas operaciones repetitivas en QGIS, Python puede proporcionarnos una solución sencilla. Antes de afrontar el trabajo deberíamos plantearnos si es posible automatizarlo de alguna manera.

En esta entrada vamos a hacer unas operaciones básicas con Python sobre una capa vectorial:

• Corrección automática de los nombres para incluir caracteres especiales.
• Adición de atributos a geometrías.
• Cálculo automático de áreas de geometrías teniendo en cuenta el Sistema de Coordenadas.

Para ello partimos de los recintos autonómicos descargados del Centro de Descargas del Instituto Geográfico Nacional (Información geográfica de referencia > Líneas límite municipales). Para cargar la capa activa en Python escribiremos la siguiente línea:

al = iface.activeLayer()

De esta forma almacenamos la capa activa en la variable al (active layer). Ahora accederemos a las entidades de la capa almacenándolas en la variable features.

features = al.getFeatures()

Supongamos que no conocemos los atributos (campos) de las entidades de la capa, mediante fields() podemos acceder a los campos de la capa. Con un bucle for los vamos sacando por cosola:

for field in al.fields():
    print(field.name())

Si quisiéramos ver los valores de cada uno de los campos para cada una de las entidades:

for feat in features:
	attrs = feat.attributes()
	print(attrs)

También podemos usar el siguiente código para mostrar el nombre del campo seguido de su valor:

for feat in features:
	for field in al.fields():
        fn = field.name()
		print(fn,": ",feat[fn])

Comprobamos que no contienen el área y que los nombres de las Comunidades Autónomas no se muestran correctamente. Para computar el área, crearemos primero ese campo en las entidades:

dp = al.dataProvider()
dp.addAttributes([QgsField("AREA", QVariant.Double)])
al.updateFields()

Pasamos la capa a dataProvider, sobre él añadimos el atributo AREA y lo definimos como de tipo decimal con QVariant.Double. Cuando acabemos forzamos mediante updateFields() los cambios en la capa. Ya solo queda calcular el área de cada una de las entidades y añadir el valor a los atributos. Para el cálculo se forzará el uso del elipsoide para la obtención de distancias (la capa cargada está en WGS84 EPSG: 4258, un sistema de coordenadas geográfico).

d = QgsDistanceArea()
d.setEllipsoid('WGS84')

Activaremos la edición de la capa. Al tiempo que iteramos sobre las entidades para obtener su área, aprovechamos para corregir sus nombres (codificados en UTF-8):

al.startEditing()
for feat in features:
    name = feat.attribute('NAMEUNIT')
    geom = feat.geometry()
    name = name.encode("windows-1252").decode("utf-8")
    feat['NAMEUNIT'] = name
    feat['AREA'] = d.convertAreaMeasurement(d.measureArea(geom), QgsUnitTypes.AreaSquareKilometers)
    al.updateFeature(feat)
    print(name)
    print("Area (km2):", d.convertAreaMeasurement(d.measureArea(geom), QgsUnitTypes.AreaSquareKilometers))
al.commitChanges()

De esta manera ya hemos añadido todos esos datos a las entidades y además corregido los nombres. Imagina que puedes hacer cálculos más complejos y sobre miles de entidades en una cuestión de minutos.

(3 votos, promedio: 5,00 de 5)

Cargando...

Formación de calidad impartida por profesionales