shapefile_area.py

# -*- coding: utf-8 -*-
"""
Created on Tue May 18 12:18:12 2010
This script takes 2 or three input parameters:

* shapefile name
* FieldName1 to aggregate on
* FieldName2 to aggregate (Optional)

It opens the shapefile, and iterates through features
aggregating the values of FieldName2 by values of FieldName1
or if FielName2 is not give, it aggregates the area of the
associated geometry by values of FieldName1.

Not tested!

@author: Jose Gomez-Dans <[email protected]>
"""

from osgeo import ogr
import sys

if len(sys.argv) <3 or len(sys.argv)>4:
    raise
try:
    ogr_dataset = ogr.Open ( sys.argv[1] )
except:
    raise IOError, "Cannot open file %s" % ( sys.argv[1] )

layer = ogr_dataset.GetLayer ( 0 )
total_area = { }
if len ( sys.argv ) == 3:
    # Only a filename and field given, calculate area based on geometry
    for feat in layer:
        lc = feat.GetFieldAsString ( feat.GetFieldName ( sys.argv[2] ) )
        geom = feat.GetGeometryRef ()
        total_area = total_area.setdefault ( lc, 0 ) + geom.GetArea() 
elif len( sys.argv ) == 4:
    #Use another column to aggregate areas
    for feat in layer:
        lc = feat.GetFieldAsString ( feat.GetFieldName ( sys.argv[2] ) )
        geom = feat.GetGeometryRef ()
        total_area = total_area.setdefault ( lc, 0 ) + \
                        feat.GetFieldAsDouble ( sys.argv[3] )

for ( lc_type, area ) in total_area.iteritems():
    print "%s\t%f\n" % ( lc_type, area )