Calcular a área total de polígonos no shapefile usando GDAL?

Eu tenho um shapefile na projeção da British National Grid:

Geometry: 3D Polygon
Feature Count: 5378
Extent: (9247.520209, 14785.170099) - (638149.173223, 1217788.569952)
Layer 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: Integer (9.0)

Posso usar GDAL / OGR para obter a área total de todos os polígonos no shapefile, em hectares?

Gostaria de saber se isso é possível com -sql, algo como:

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

Mas tentando isso eu entendo ERROR 1: Undefined function 'ST_Area' used..

Acho que poderia importar o Shapefile para o QGIS, adicionar um atributo de área a cada polígono e depois somar, mas prefiro usar uma ferramenta de linha de comando, se possível.

gdal ogr area ogrinfo Richard
fonte

Respostas:

Há um campo especial no OGR SQL chamado OGR_GEOM_AREA que retorna a área da geometria do recurso:

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

onde a TOTAL_AREAunidade de medida depende da camada SRS (leia os comentários abaixo).

Antonio Falciano
fonte

Surpreendente! Isso me dá SUM_OGR_GEOM_AREA (Real) = 4459037129.50955. Isso está em hectares ou em alguma outra unidade? E importa em que projeção está o meu shapefile de origem?

Richard

A unidade de medida deriva de como suas geometrias são projetadas e, portanto, são metros quadrados.

Antonio Falciano 21/03

mais provável em m ^ 2, divisão por 10000 para converter a ha

DMCI

Obrigado! Acabei de encontrar os documentos: gdal.org/classOGRSurface.html#a3b2c3125ec8c0b3a986e43cd1056f9e4 Retorna a área do recurso em unidades quadradas do sistema de referência espacial em uso . Desculpe por ser um iniciante, mas, caso esteja usando um shapefile em um SRS diferente, como descubro quais são as unidades quadradas de um SRS específico?

Richard

Olhando em sua camada SRS WKT (ou seja, o arquivo de projeção), você encontrará UNIT ["Meter", 1]]

Antonio Falciano

Sim, é possível, mas você precisa usar o dialeto OGR SQLite da seguinte maneira:

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

Além disso, verifique se esse myshapefileé o nome da camada myshapefile.shp. Você pode fazer isso da seguinte maneira:

ogrinfo myshapefile.shp

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

dmci
fonte

Obrigado! Na verdade eu vejo no such column: geometry. Como descubro como é chamada a coluna de geometria? Esta é a saída ogrinfo -al -fid 1:

OGRFeature(mylayer):1   cat (Integer) = 2   POLYGON ((463267.036276041297242 1216886.583904854720458 0,463267.693611663184129 1216956.525473011657596 0,463405.369117364054546 1216820.109560555079952 0,463404.712737055611797 1216750.1665881925728170,463267.036276041297242 1216886.583904854720458 0,463267.036276041297242 1216886.583904854720458 0))

Richard

Do ogrinfo -so -alBut @dmci, eu não entendo nada na sua resposta: pois os shapefiles da GDAL sempre têm uma camada e seu nome é o nome base do shapefile. Como você pode obter o nome "mytable" em vez de "myshapefile"?

precisa saber é o seguinte

Obrigado! A saída desse comando está na pergunta original.

Richard

@ user30184 - concordo totalmente - mas eu estava replicando as convenções de nomenclatura que o OP havia usado na pergunta. Para

manter a

Ok, parece ser específico do driver. Alguns drivers relatam isso Geometry Column = GEOMETRYcom o ogrinfo -so -al, mas o driver shapefile não. Com os shapefiles, suponho que o nome seja "OGR_GEOMETRY" para o dialeto OGR SQL (consulte os campos especiais em gdal.org/ogr_sql.html ) e "geometria" para o dialeto SQLite.

user30184

Usando o QGIS, você pode executar este código simples para imprimir a área total do shapefile (suponho que você esteja avaliando a área em um sistema de referência projetado):

from qgis.core import *

filepath = 'C:/Users/path_to_the_shapefile/shapefile.shp'
layer = QgsVectorLayer(filepath, 'layer' , 'ogr')

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

print area # total area in square meters

print area/10000 # total area in hectares

mgri
fonte