11

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.

15

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).

  • 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 '17 at 11:05
  • 1
    The unit of measure derives from how your geometries are projected, so they are square meters. – Antonio Falciano Mar 21 '17 at 11:08
  • 1
    most likely in m^2, divide by 10000 to convert to hectares – dmci Mar 21 '17 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 '17 at 11:10
  • 1
    Looking in your layer SRS WKT (i.e. the projection file), you will find UNIT["Meter",1]] – Antonio Falciano Mar 21 '17 at 11:13
6

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)
  • 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 '17 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 '17 at 11:04
  • Thanks! The output from that command is in the original question. – Richard Mar 21 '17 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 '17 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 '17 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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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