Calculer la superficie totale des polygones dans le fichier de formes à l'aide de GDAL?

11

J'ai un fichier de formes dans la projection 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)

Puis-je utiliser GDAL / OGR pour obtenir la superficie totale de tous les polygones du fichier de formes, en hectares?

Je me demande si c'est possible avec -sqlquelque chose comme:

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

Mais j'essaye ERROR 1: Undefined function 'ST_Area' used..

Je suppose que je pourrais importer le Shapefile dans QGIS, ajouter un attribut de zone à chaque polygone, puis le résumer, mais je préférerais de loin utiliser un outil de ligne de commande si possible.

gdal ogr area ogrinfo Richard
la source

Réponses:

15

Il y a un champ spécial dans OGR SQL appelé OGR_GEOM_AREA qui renvoie la zone de la géométrie de l' entité :

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

où l' TOTAL_AREAunité de mesure dépend de la couche SRS (lire les commentaires ci-dessous).

Antonio Falciano
la source
Incroyable! Cela me donne SUM_OGR_GEOM_AREA (Real) = 4459037129.50955. Est-ce en hectares ou dans une autre unité? Et quelle importance a la projection de mon fichier de formes source?
Richard
1
L'unité de mesure découle de la façon dont vos géométries sont projetées, elles sont donc en mètres carrés.
Antonio Falciano
1
très probablement en m ^ 2, divisez par 10000 pour convertir en hectares
dmci
Merci! Je viens de trouver les documents: gdal.org/classOGRSurface.html#a3b2c3125ec8c0b3a986e43cd1056f9e4 Renvoie la zone de l' entité en unités carrées du système de référence spatiale utilisé . Désolé d'être un débutant total, mais au cas où j'utiliserais un fichier de formes dans un SRS différent, comment savoir quelles sont les unités carrées d'un SRS particulier?
Richard
1
En regardant dans votre couche SRS WKT (ie le fichier de projection), vous trouverez UNIT ["Meter", 1]]
Antonio Falciano
6

Oui, c'est possible, mais vous devez utiliser le dialecte OGR SQLite comme suit:

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

Assurez-vous également qu'il myshapefiles'agit du nom de la couche myshapefile.shp. Vous pouvez le faire comme suit:

ogrinfo myshapefile.shp

INFO: Open of `myshapefile.shp`
    using driver `ESRI Shapefile' successful.
1: myshapefile (Polygon)
dmci
la source
Merci! En fait, je vois no such column: geometry. Comment savoir comment s'appelle la colonne de géométrie? Ceci est la sortie de 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
1
Ne ogrinfo -so -alfaites pas @dmci, je ne comprends pas un bit dans votre réponse: pour les fichiers de formes GDAL ont toujours une couche et son nom est le nom de base du fichier de formes. Comment obtenir le nom "mytable" au lieu de "myshapefile"?
user30184
Merci! La sortie de cette commande se trouve dans la question d'origine.
Richard
@ user30184 - tout à fait d'accord - mais je reproduisais les conventions de dénomination que l'OP avait utilisées dans leur question. Par
souci de
2
D'accord, cela semble être spécifique au pilote. Certains pilotes le signalent comme Geometry Column = GEOMETRYavec ogrinfo -so -al mais pas le pilote shapefile. Avec les fichiers de formes, je suppose que le nom est "OGR_GEOMETRY" pour le dialecte SQL OGR (voir les champs spéciaux dans gdal.org/ogr_sql.html ) et "geometry" pour le dialecte SQLite.
user30184
4

En utilisant QGIS, vous pouvez exécuter ce code simple pour imprimer la zone totale du fichier de formes (je suppose que vous évaluez la zone dans un système de référence projeté):

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
la source