Oblicz całkowity obszar wielokątów w pliku kształtu za pomocą GDAL?

11

Mam plik kształtu w projekcji 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)

Czy mogę użyć GDAL / OGR, aby uzyskać całkowitą powierzchnię wszystkich wielokątów w pliku kształtu, w hektarach?

Zastanawiam się, czy jest to możliwe przy pomocy -sqlczegoś takiego:

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

Ale próbuję to ERROR 1: Undefined function 'ST_Area' used..

Myślę, że mógłbym zaimportować plik Shapefile do QGIS, dodać atrybut obszaru do każdego wielokąta, a następnie zsumować go, ale wolałbym raczej użyć narzędzia wiersza poleceń, jeśli to możliwe.

gdal ogr area ogrinfo Richard
źródło

Odpowiedzi:

15

W OGR SQL znajduje się specjalne pole o nazwie OGR_GEOM_AREA, które zwraca obszar geometrii elementu:

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

gdzie TOTAL_AREAjednostka miary zależy od warstwy SRS (przeczytaj komentarze poniżej).

Antonio Falciano
źródło
Niesamowity! To mi daje SUM_OGR_GEOM_AREA (Real) = 4459037129.50955. Czy to w hektarach, czy w innej jednostce? I czy ma to znaczenie, w jakiej projekcji znajduje się mój plik kształtu?
Richard
1
Jednostka miary wywodzi się z rzutu twoich geometrii, więc są to metry kwadratowe.
Antonio Falciano
1
najprawdopodobniej wm ^ 2, podziel przez 10000, aby przeliczyć na hektary
dmci
Dzięki! Właśnie znalazłem dokumentację: gdal.org/classOGRSurface.html#a3b2c3125ec8c0b3a986e43cd1056f9e4 Zwraca obszar funkcji w jednostkach kwadratowych używanego systemu odniesienia przestrzennego . Przepraszam, że jestem całkowicie początkującym, ale na wypadek, gdyby korzystałem z pliku shapefile w innym SRS, jak mogę dowiedzieć się, jakie są jednostki kwadratowe w konkretnym SRS?
Richard
1
Przyglądając się swojej warstwie SRS WKT (tj. Plik projekcji), znajdziesz UNIT [„Meter”, 1]]
Antonio Falciano
6

Tak, jest to możliwe, ale musisz użyć dialektu OGR SQLite w następujący sposób:

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

Upewnij się również, że myshapefilejest to nazwa warstwy w myshapefile.shp. Możesz to zrobić w następujący sposób:

ogrinfo myshapefile.shp

INFO: Open of `myshapefile.shp`
    using driver `ESRI Shapefile' successful.
1: myshapefile (Polygon)
dmci
źródło
Dzięki! Właściwie to rozumiem no such column: geometry. Jak dowiedzieć się, jak nazywa się kolumna geometrii? Jest to wyjście 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, nie rozumiem ani odrobiny w twojej odpowiedzi: ponieważ pliki kształtów GDAL mają zawsze jedną warstwę, a jej nazwa to nazwa_pliku pliku kształtu. Jak uzyskać nazwę „mytable” zamiast „myshapefile”?
user30184
Dzięki! Dane wyjściowe tego polecenia znajdują się w pierwotnym pytaniu.
Richard
@ user30184 - całkowicie się zgadzam - ale replikowałem konwencje nazewnictwa, których użył OP w swoim pytaniu. Dla spójności zaktualizowałem własną odpowiedź
dmci
2
Ok, wydaje się, że jest specyficzny dla sterownika. Niektórzy kierowcy zgłaszają to tak, jak w Geometry Column = GEOMETRYprzypadku ogrinfo -so -al, ale sterownik shapefile nie. W przypadku plików kształtów przypuszczam, że nazwa to „OGR_GEOMETRY” dla dialektu SQL OGR (patrz specjalne pola w gdal.org/ogr_sql.html ) i „geometria” dla dialektu SQLite.
user30184
4

Korzystając z QGIS, możesz uruchomić ten prosty kod do wydrukowania całkowitego obszaru pliku shapefile (zakładam, że oceniasz obszar w rzutowanym systemie odniesienia):

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
źródło