Hitung luas total poligon di shapefile menggunakan GDAL?

11

Saya memiliki shapefile dalam proyeksi 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)

Bisakah saya menggunakan GDAL / OGR untuk mendapatkan luas total semua poligon di shapefile, dalam hektar?

Saya ingin tahu apakah ini mungkin terjadi -sql, seperti:

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

Tetapi mencoba itu saya dapatkan ERROR 1: Undefined function 'ST_Area' used..

Saya kira saya bisa mengimpor Shapefile ke QGIS, menambahkan atribut area ke setiap poligon, dan kemudian menjumlahkannya, tetapi saya lebih suka menggunakan alat baris perintah jika memungkinkan.

gdal ogr area ogrinfo Richard
sumber

Jawaban:

15

Ada bidang khusus dalam OGR SQL yang disebut OGR_GEOM_AREA yang mengembalikan bidang geometri fitur:

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

di mana TOTAL_AREAsatuan ukuran tergantung pada lapisan SRS (baca komentar di bawah).

Antonio Falciano
sumber
Luar biasa! Ini memberiku SUM_OGR_GEOM_AREA (Real) = 4459037129.50955. Apakah ini dalam hektar atau unit lain? Dan apakah penting proyeksi apa yang ada di formfile sumber saya?
Richard
1
Satuan ukuran berasal dari bagaimana geometri Anda diproyeksikan, jadi itu adalah meter persegi.
Antonio Falciano
1
kemungkinan besar dalam m ^ 2, dibagi dengan 10.000 untuk dikonversi menjadi hektar
dmci
Terima kasih! Saya baru saja menemukan dokumen: gdal.org/classOGRSurface.html#a3b2c3125ec8c0b3a986e43cd1056f9e4 Mengembalikan area fitur dalam satuan kuadrat dari sistem referensi spasial yang digunakan . Maaf menjadi pemula total, tetapi untuk berjaga-jaga jika saya menggunakan shapefile di SRS yang berbeda, bagaimana cara mengetahui apa satuan kuadrat dari SRS tertentu?
Richard
1
Melihat di lapisan Anda SRS WKT (yaitu file proyeksi), Anda akan menemukan UNIT ["Meter", 1]]
Antonio Falciano
6

Ya, itu mungkin, tetapi Anda perlu menggunakan dialek OGR SQLite sebagai berikut:

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

Juga, pastikan itu myshapefileadalah nama layer di myshapefile.shp. Anda dapat melakukan ini sebagai berikut:

ogrinfo myshapefile.shp

INFO: Open of `myshapefile.shp`
    using driver `ESRI Shapefile' successful.
1: myshapefile (Polygon)
dmci
sumber
Terima kasih! Sebenarnya saya mengerti no such column: geometry. Bagaimana cara mengetahui apa yang disebut kolom geometri? Ini adalah output dari 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 -alBut @dmci, saya tidak mengerti sedikit pun dalam jawaban Anda: untuk GDAL shapefile selalu memiliki satu layer dan namanya adalah nama dasar dari shapefile. Bagaimana Anda bisa mendapatkan nama "mytable" daripada "myshapefile"?
user30184
Terima kasih! Output dari perintah itu ada di pertanyaan awal.
Richard
@ user30184 - sepenuhnya setuju - tapi saya mereplikasi konvensi penamaan yang telah digunakan OP dalam pertanyaan mereka. Untuk konsistensi, saya telah memperbarui jawaban saya sendiri
dmci
2
Ok, sepertinya driver khusus. Beberapa driver melaporkannya seperti Geometry Column = GEOMETRYdengan ogrinfo -jadi -al tetapi driver shapefile tidak. Dengan shapefile saya kira namanya adalah "OGR_GEOMETRY" untuk dialek OGR SQL (lihat bidang khusus di gdal.org/ogr_sql.html ) dan "geometri" untuk dialek SQLite.
user30184
4

Dengan menggunakan QGIS, Anda dapat menjalankan kode sederhana ini untuk mencetak total area shapefile (saya berasumsi Anda sedang mengevaluasi area dalam sistem referensi yang diproyeksikan):

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
sumber