Calcular el área total de los polígonos en el shapefile usando GDAL.

Question

Tengo un shapefile en proyección British National Grid:

Geometría: Polígono 3D
Recuento de entidades: 5378
Extensión: (9247.520209, 14785.170099) - (638149.173223, 1217788.569952)
Capa SRS WKT:
PROJCS["British_National_Grid",
    GEOGCS["GCS_airy",
        DATUM["OSGB_1936",
            SPHEROID["Airy_1830",6377563.396,299.3249646]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",49],
    PARAMETER["central_meridian",-2],
    PARAMETER["scale_factor",0.9996012717],
    PARAMETER["false_easting",400000],
    PARAMETER["false_northing",-100000],
    UNIT["Meter",1]]
cat: Entero (9.0)

¿Puedo usar GDAL/OGR para obtener el área total de todos los polígonos en el shapefile, en hectáreas?

Me pregunto si esto es posible con -sql, algo así como:

ogrinfo -sql "SELECT SUM(ST_Area(geom::geography)) FROM mytable" myshapefile.shp

Pero al intentarlo obtengo ERROR 1: Undefined function 'ST_Area' used..

Supongo que podría importar el Shapefile en QGIS, agregar un atributo de área a cada polígono y luego sumarlo, pero preferiría usar una herramienta de línea de comandos si es posible.

Answer 1

Existe un campo especial en OGR SQL llamado OGR_GEOM_AREA que devuelve el área de la geometría de la entidad:

ogrinfo -sql "SELECT SUM(OGR_GEOM_AREA) AS TOTAL_AREA FROM myshapefile" myshapefile.shp

donde la unidad de medida de TOTAL_AREA depende del SRS de la capa (lee los comentarios abajo).

Answer 2

Sí, es posible, pero necesitas usar el dialecto SQLite de OGR como sigue:

ogrinfo -dialect SQLite -sql 'SELECT SUM(ST_Area(geometry))/10000 FROM myshapefile' myshapefile.shp

También asegúrate de que myshapefile sea el nombre de la capa en myshapefile.shp. Puedes hacer esto de la siguiente manera:

ogrinfo myshapefile.shp

INFO: Open of `myshapefile.shp`
    using driver `ESRI Shapefile' successful.
1: myshapefile (Polygon)

Answer 3

Usando QGIS, puedes ejecutar este código simple para imprimir el área total del archivo de formas (asumo que estás evaluando el área en un sistema de referencia proyectado):

from qgis.core import *

filepath = 'C:/Usuarios/ruta_al_archivo_de_formas/archivo_de_formas.shp'
capa = QgsVectorLayer(filepath, 'capa' , 'ogr')

area = 0
for feat in capa.getFeatures():
    area += feat.geometry().area()

print area # área total en metros cuadrados

print area/10000 # área total en hectáreas