OGR and the area of a polygon

· Read in about 2 min · (278 Words)

I was writing a little utility to process a large number of polygons in a bunch of shape files, and wanted to have the area of each polygon. I thought I will do it using Python and OGR library. All my files are in WGS84 Latitude Longitude projection.

For experiment I tried the following code:

from osgeo import ogr

driver = ogr.GetDriverByName('ESRI Shapefile')
hinputfile = driver.Open(r'c:\GeoData\AZ.shp', 0)
in_layer = hinputfile.GetLayer(0)
feature = in_layer.GetFeature(0)
geom = feature.GetGeometryRef()
print(geom.GetArea())

This would give the area in square degrees, which I definitely not want. This would not be a good indicator, since the actual area would get smaller and smaller for the same value the bigger the latitude gets. I wanted to get the area in m² or km². This can be easily achieved by using a projection with units as meters. For example the «google» projection (EPSG:900913 - gOOglE get it?) which is now EPSG:3857 which is in meters. Here is the WKT for it.

PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["Popular Visualisation CRS",
        DATUM["Popular_Visualisation_Datum",
            SPHEROID["Popular Visualisation Sphere",6378137,0,
                AUTHORITY["EPSG","7059"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6055"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.01745329251994328,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4055"]],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    AUTHORITY["EPSG","3785"],
    AXIS["X",EAST],
    AXIS["Y",NORTH]]

So the code to calculate the area of each polygon in a shapefile in sq km would look something like the following:

from osgeo import ogr
from osgeo import osr

driver = ogr.GetDriverByName('ESRI Shapefile')

hinputfile = driver.Open(r'c:\GeoData\AZ.shp', 0)
in_layer = hinputfile.GetLayer(0)

# Prepare a transformation between projections
src_srs = in_layer.GetSpatialRef()
tgt_srs = osr.SpatialReference()
tgt_srs.ImportFromEPSG(3857)
transform = osr.CoordinateTransformation(src_srs, tgt_srs)

for feature in in_layer:
    geom = feature.GetGeometryRef()

    geom2 = geom.Clone()
    geom2.Transform(transform)

    area_in_sq_m = geom2.GetArea()
    area_in_sq_km = area_in_sq_m / 1000000
    print('Area in sq km: ', area_in_sq_km)

hinputfile.Destroy()

Hope this helps.