Ejercicio con funciones en PyQGIS: determinación de áreas

Como ya se ha dicho en otra oportunidad en relación al código PyQGIS que se encuentra en diferentes posts de este Blog, la mayoría de ellos se limita a probar en la Python Console algunos métodos de clase con diferentes procedimientos que realizan tareas útiles empleando archivos vectoriales o ráster. Ello es así porque existen más de 1200 clases en este entorno de programación y se trata de encontrar procedimientos que permitan introducir las clases más representativas y la forma de explorar sus métodos utilizando los recursos y la documentación que pone a la disposición el propio sistema. Es por ello que se ha insistido en el uso de la función ‘getPat’ y el plugin ‘Get Pattern’.

Por otra parte, es preciso añadir que gran parte de esa información está particularizada y su utilidad se vería incrementada si se generaliza (para evitar errores en tiempo de ejecución) y pone en forma de función. Ello permitiría la adaptación y reutilización del código; sobre todo si pensamos en términos de plugins futuros.

Supongamos que estamos interesados en crear una aplicación que a partir del rasgo seleccionado de un vectorial tipo polígono, con una proyección en grados, se determine el área correspondiente. Ese problema está contemplado en este Blog en el post:

Áreas de polígonos en sistemas proyectados en coordenadas geográficas con PyQGIS

Sin embargo, el código está particularizado a un sólo rasgo y me interesaría generalizarlo en forma de función para que tome en cuenta lo siguiente:

  • 1) Comprobar si existe capa activa e interrumpir la ejecución en caso negativo.
  • 2) Comprobar el crs de la capa e intentar hacer la determinación del área sólo si éste está en coordenadas geográficas.
  • 3) Comprobar si existen rasgos seleccionados y en caso afirmativo determinar el área (en km2) para todos ellos.
  • 4) Comprobar el tipo de geometría de un rasgo e interrumpir la ejecución si éste no es del tipo polígono.

En el caso del numeral 3, sabemos que un objeto de la clase QgsVectorLayer se convierte en un objeto de la clase QgsCoordinateReferenceSystem mediante el método ‘crs’ de la primera. En esta última, el método ‘geographicFlag’ permite entonces hacer la comprobación produciendo como resultado ‘False’ o ‘True’.

Por otra parte, con relación al numeral 4, tenemos que algunas comprobaciones se realizan frente a los tipos enumerados que representan diferentes objetos en las clases. Estos tipos enumerados se encuentran concentrados en la clase QGis y permiten que el código sea más legible porque la comprobación se realiza escribiendo el texto correspondiente al objeto de tipo enumerado y no un número. Además, nos ahorra el trabajo de crear un diccionario. Por ejemplo, sabemos que una geometría de tipo polígono tiene por valor númerico el 2. Si se quiere averiguar por código el tipo de geometría de un objeto QgsVectorLayer entonces: layer.feature().geometry().type(). Si es de tipo polígono el valor devuelto es 2. La sintaxis que equivale a ese valor es QGis.Polygon y es la que se debería usar en la comprobación; no el 2.

Tomando en cuenta lo anterior, si queremos resumir el código del post considerado arriba se podría proponer lo siguiente para la función:

from qgis.gui import QgsMessageBar

def area_geo():

    layer = iface.activeLayer()
    
    if layer is None or (layer.type() == QgsMapLayer.RasterLayer):
        iface.messageBar().pushMessage("Warning:",
                                       u"There is not Vector Active Layer",
                                       QgsMessageBar.WARNING, 5) 
        return
    
    crs = layer.crs()
    
    feats = layer.selectedFeatures()

    if len(feats) == 0:
        iface.messageBar().pushMessage("Warning:",
                                       u"There is not selected features",
                                       QgsMessageBar.WARNING, 5) 
        
    else:
    
        if crs.geographicFlag() == True:
            
            objectArea = QgsDistanceArea()    #creando objeto
            objectArea.setEllipsoid('WGS84')  #estableciendo elipsoide
            objectArea.setEllipsoidalMode(True)  #estableciendo modo elipsoidal

            for feature in feats:
                if feature.geometry().type() != QGis.Polygon:
                    iface.messageBar().pushMessage("Warning:",
                                                    u"Geometry is not a Polygon",
                                                    QgsMessageBar.WARNING, 5)
                    return
                qgsPointsList = feature.geometry().asPolygon()
                area = objectArea.measurePolygon(qgsPointsList[0])   #determinando area
                print "area = %.2f km2" % (area/1e6)
        
        else:
            iface.messageBar().pushMessage("Warning:",
                                            u"CRS is not in geographic coordinates",
                                            QgsMessageBar.WARNING, 5)

Al ejecutarla sin ninguna capa activa en el map Canvas:

function1

Si cargamos el vectorial world_borders (EPSG: 4326); pero sin seleccionar ningún rasgo para determinar el área:

function3

Finalmente, en esta última situación se han seleccionado los rasgos continentales (excluyendo los insulares) de Italia, Portugal, Francia y España. Los valores de las áreas que refleja la Python Console son, respectivamente, las de los paises mencionados con anterioridad.

function2

Esta entrada fue publicada en PyQGIS, QGIS, Software Libre. Guarda el enlace permanente.

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s