14

I have a shapefile in British National Grid projection:

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)

Can I use GDAL/OGR to get the total area of all the polygons in the shapefile, in hectares?

I'm wondering if this is possible with -sql, something like:

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

But trying that I get ERROR 1: Undefined function 'ST_Area' used..

I guess I could import the Shapefile into QGIS, add an area attribute to each polygon, and then sum it, but I would much rather use a command line tool if possible.

Improve this question

3 Answers 3

Reset to default
19

There's a special field in OGR SQL called OGR_GEOM_AREA which returns the area of the feature's geometry:

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

where TOTAL_AREA unit of measure depends by the layer SRS (read the comments below).

Improve this answer
5
  • 1
    Amazing! This gives me SUM_OGR_GEOM_AREA (Real) = 4459037129.50955. Is this in hectares or some other unit? And does it matter what projection my source shapefile is in?
    – Richard
    Mar 21, 2017 at 11:05
  • 2
    The unit of measure derives from how your geometries are projected, so they are square meters.
    – Antonio Falciano
    Mar 21, 2017 at 11:08
  • 2
    most likely in m^2, divide by 10000 to convert to hectares
    – dmci
    Mar 21, 2017 at 11:08
  • Thanks! I just found the docs: gdal.org/classOGRSurface.html#a3b2c3125ec8c0b3a986e43cd1056f9e4 Returns the area of the feature in square units of the spatial reference system in use. Sorry to be a total beginner, but just in case I'm using a shapefile in a different SRS, how do I find out what the square units of a particular SRS are?
    – Richard
    Mar 21, 2017 at 11:10
  • 2
    Looking in your layer SRS WKT (i.e. the projection file), you will find UNIT["Meter",1]]
    – Antonio Falciano
    Mar 21, 2017 at 11:13
11

Yes, it is possible, but you need to use the OGR SQLite dialect as follows:

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

Also, ensure that myshapefile is the layer name in myshapefile.shp. You can do this as follows:

ogrinfo myshapefile.shp

INFO: Open of `myshapefile.shp`
    using driver `ESRI Shapefile' successful.
1: myshapefile (Polygon)
Improve this answer
6
  • Thanks! Actually I see no such column: geometry. How do I find out what the geometry column is called? This is the output of 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
    Mar 21, 2017 at 11:03
  • 1
    Do ogrinfo -so -al But @dmci, I do not understand one bit in your answer: for GDAL shapefiles have always one layer and its name is the basename of the shapefile. How can you get name "mytable" instead of "myshapefile"?
    – user30184
    Mar 21, 2017 at 11:04
  • Thanks! The output from that command is in the original question.
    – Richard
    Mar 21, 2017 at 11:05
  • @user30184 - totally agree - but I was replicating the naming conventions that the OP had used in their question. For consistency, I've updated my own answer
    – dmci
    Mar 21, 2017 at 20:58
  • 2
    Ok, it seems to be driver specific. Some drivers report it like Geometry Column = GEOMETRY with ogrinfo -so -al but shapefile driver does not. With shapefiles I suppose that the name is "OGR_GEOMETRY" for OGR SQL dialect (see special fields in gdal.org/ogr_sql.html) and "geometry" for SQLite dialect.
    – user30184
    Mar 22, 2017 at 9:09
4

Using QGIS, you may run this simple code for printing the total area of the shapefile (I assume you are evaluating the area in a projected reference system):

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
Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.