Gesamtfläche von Polygonen im Shapefile mit GDAL berechnen?

11

Ich habe ein Shapefile in der British National Grid-Projektion:

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)

Kann ich GDAL / OGR verwenden, um die Gesamtfläche aller Polygone im Shapefile in Hektar zu ermitteln?

Ich frage mich, ob dies möglich ist mit -sql:

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

Aber das versuche ich zu bekommen ERROR 1: Undefined function 'ST_Area' used..

Ich denke, ich könnte das Shapefile in QGIS importieren, jedem Polygon ein Bereichsattribut hinzufügen und es dann summieren, aber ich würde viel lieber ein Befehlszeilen-Tool verwenden, wenn dies möglich ist.

gdal ogr area ogrinfo Richard
quelle

Antworten:

15

In OGR SQL gibt es ein spezielles Feld namens OGR_GEOM_AREA, das den Bereich der Geometrie des Features zurückgibt:

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

wobei die TOTAL_AREAMaßeinheit von der Schicht SRS abhängt (lesen Sie die Kommentare unten).

Antonio Falciano
quelle
Tolle! Das gibt mir SUM_OGR_GEOM_AREA (Real) = 4459037129.50955. Ist das in Hektar oder einer anderen Einheit? Und spielt es eine Rolle, in welcher Projektion sich mein Quell-Shapefile befindet?
Richard
1
Die Maßeinheit ergibt sich aus der Projektion Ihrer Geometrien, also Quadratmeter.
Antonio Falciano
1
höchstwahrscheinlich in m ^ 2 durch 10000 dividieren, um in Hektar
umzurechnen
Vielen Dank! Ich habe gerade die Dokumente gefunden: gdal.org/classOGRSurface.html#a3b2c3125ec8c0b3a986e43cd1056f9e4 Gibt den Bereich des Features in quadratischen Einheiten des verwendeten Raumbezugssystems zurück . Es tut mir leid, ein absoluter Anfänger zu sein, aber für den Fall, dass ich ein Shapefile in einem anderen SRS verwende, wie finde ich heraus, wie groß die quadratischen Einheiten eines bestimmten SRS sind?
Richard
1
Wenn Sie in Ihre Ebene SRS WKT (dh die Projektionsdatei) schauen, finden Sie UNIT ["Meter", 1]]
Antonio Falciano
6

Ja, es ist möglich, aber Sie müssen den OGR SQLite-Dialekt wie folgt verwenden:

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

Stellen Sie außerdem sicher, dass dies myshapefileder Layername in ist myshapefile.shp. Sie können dies wie folgt tun:

ogrinfo myshapefile.shp

INFO: Open of `myshapefile.shp`
    using driver `ESRI Shapefile' successful.
1: myshapefile (Polygon)
dmci
quelle
Vielen Dank! Eigentlich verstehe ich no such column: geometry. Wie finde ich heraus, wie die Geometriespalte heißt? Dies ist die Ausgabe von 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
Do ogrinfo -so -alAber @dmci, ich verstehe kein bisschen in Ihrer Antwort: Für GDAL-Shapefiles gibt es immer eine Ebene und der Name ist der Basisname des Shapefiles. Wie kann man den Namen "mytable" anstelle von "myshapefile" bekommen?
user30184
Vielen Dank! Die Ausgabe dieses Befehls befindet sich in der ursprünglichen Frage.
Richard
@ user30184 - stimme vollkommen zu - aber ich habe die Namenskonventionen repliziert, die das OP in seiner Frage verwendet hatte. Aus
Gründen der
2
Ok, es scheint fahrerspezifisch zu sein. Einige Treiber melden dies wie Geometry Column = GEOMETRYbei ogrinfo -so -al, Shapefile-Treiber jedoch nicht. Bei Shapefiles ist der Name "OGR_GEOMETRY" für den OGR-SQL-Dialekt (siehe spezielle Felder in gdal.org/ogr_sql.html ) und "Geometrie" für den SQLite-Dialekt.
user30184
4

Mit QGIS können Sie diesen einfachen Code zum Drucken der Gesamtfläche des Shapefiles ausführen (ich gehe davon aus, dass Sie die Fläche in einem projizierten Referenzsystem auswerten):

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
quelle