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

Normalmente, en sistemas proyectados con coordenadas geográficas se hace mención de que no se pueden determinar elementos de geometría porque la proyección no está en metros. Sin embargo, mediante PyQGIS se puede a través de la clase QgsDistanceArea. Estos cálculos se basan en el elipsoide y las coordenadas de los puntos tienen que estar en coordenadas geográficas. En esta oportunidad vamos a calcular el área del territorio continental de España a partir de un mapa de circulación libre en Internet denominado world_borders.shp.

Cargando el shapefile a QGIS, hacemos click en España teniendo habilitada la herramienta de “Select Single Feature” (Attributes Tool Bar) y luego en “Zoom To Selection” (Map Navigation Tool Bar). El resultado luce de la siguiente manera:

area

En la Python Console hacemos un dir(QgsDistanceArea) para explorar sus métodos:

>>>dir(QgsDistanceArea)
['__class__', '__delattr__', '__dict__', '__doc__', '__format__', 
'__getattribute__', '__hash__', '__init__', '__module__', 
'__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', 
'__sizeof__', '__str__', '__subclasshook__', '__weakref__', 
'bearing', 'convertMeasurement', 'ellipsoid', 'ellipsoidInverseFlattening', 
'ellipsoidSemiMajor', 'ellipsoidSemiMinor', 'ellipsoidalEnabled', 
'geographic', 'measure', 'measureLine', 'measurePerimeter', 
'measurePolygon', 'setEllipsoid', 'setEllipsoidalMode', 'setSourceAuthId', 
'setSourceCrs', 'sourceCrs', 'textUnit']

El método que determina el área es ‘measurePolygon’ que acepta como argumentos una lista de QgsPoint en coordenadas geográficas. Para obtener los QgsPoint del polígono correspondiente a España tenemos que obtener primero su ‘feature’, es decir, un objeto de la clase QgsFeature. Esto se logra con el método ‘selectedFeatures’ de QgsVectorLayer. Luego se convierte en un objeto de la clase QgsGeometry para utilizar su método ‘asPolygon’ que nos proporciona la lista de QgsPoint. El procedimiento sería:

>>>wb=iface.activeLayer()
>>>feat = wb.selectedFeatures()
>>>feat
[<qgis.core.QgsFeature object at 0x20B4B6F0>]
>>>geom = feat[0].geometry().asPolygon()
>>>len(geom[0])  #length de geom
1143  #QgsPoint
>>>geom[0][0:10]  #slice de geom con 11 elementos
[(-0.761016,37.846), (-0.801389,37.8014), (-0.841944,37.7489), 
(-0.858333,37.7217), (-0.86,37.7167), (-0.8575,37.7119), 
(-0.814167,37.6647), (-0.804722,37.6572), (-0.786945,37.6478), 
(-0.740833,37.6303)]
>>>geom[0][-1]
(-0.761016,37.846)   #ultimo punto que coincide con el primero

Una vez obtenida la lista de QgsPoint ahora resta crear el objeto (clase QgsDistanceArea), establecer el elipsoide (‘setEllipsoid’) y el modo elipsoidal (‘setEllipsoidalMode’) antes de usar ‘measurePolygon’. Por tanto:

>>>area = QgsDistanceArea()    #creando objeto
>>>area.setEllipsoid('WGS84')  #estableciendo elipsoide
True
>>>area.setEllipsoidalMode(True)  #estableciendo modo elipsoidal
>>>area.measurePolygon(geom[0])   #determinando area
494564854638.19434  #metros cuadrados

El area de la porción continental de España es de 494.564,85 km2.

Esta entrada fue publicada en Código Python, PyQGIS, SIG, Software Libre. Guarda el enlace permanente.

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

  1. geoariel dijo:

    Cuando introduzco el siguiente código:

    area.setEllipsoidalMode(True)

    Me entrega el siguiente error:

    Traceback (most recent call last):
    File “”, line 1, in
    AttributeError: ‘QgsDistanceArea’ object has no attribute ‘setEllipsoidal’

  2. Gracias por el comentario pero ya la Python Console te había contestado:

    AttributeError: ‘QgsDistanceArea’ object has no attribute ‘setEllipsoidal’

    Obviaste escribir el resto del método. Te faltó ‘Mode’. Es ‘setEllipsoidalMode’.

    • geoariel dijo:

      Como se ve en mi comentario, eso fue lo que puse (area.setEllipsoidalMode(True)). Sin embargo, me generó tal error.

      • Si esto:

        AttributeError: ‘QgsDistanceArea’ object has no attribute ‘setEllipsoidal’

        fue lo que te devolvió la Python Console significa que no lo hiciste como debía. La Python Console “tiene razón”. Ese método no existe en esa clase.

  3. Pingback: Áreas de polígonos en sistemas pr...

  4. Pingback: Ejercicio con funciones en PyQGIS: determinación de áreas | El Blog de José Guerrero

  5. Pingback: Determinación de longitudes en sistemas proyectados en coordenadas geográficas con PyQGIS | El Blog de José Guerrero

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