gdal - Calculate area of polygons using OGR in python script

I have a shapefile of many thousands of polygons. I'm trying to append an area field to the attribute table and calculate the area of each polygon (in sq. meters)

driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(path_to_shp_data, 1)
layer = dataSource.GetLayer()
new_field = ogr.FieldDefn("Area", ogr.OFTReal)
new_field.SetWidth(32)
layer.CreateField(new_field)


for feature in layer:
    geom = feature.GetGeometryRef()
    area = geom.GetArea() 
    feature.SetField("Area", area)
    layer.SetFeature(feature)

dataSource = None

The code runs successfully but the resulting "Area" field contains all 0's.

Projection info is as follows:

PROJCS["WGS_1984_UTM_Zone_48N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",105],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]

Answer

I ran your script (slightly modified) at the Python Console of QGIS:

from osgeo import ogr

vlayer = iface.activeLayer()

provider = vlayer.dataProvider()


path = provider.dataSourceUri()

tmp = path.split("|")

path_to_shp_data = tmp[0]

driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(path_to_shp_data, 1)
layer = dataSource.GetLayer()
new_field = ogr.FieldDefn("Area", ogr.OFTReal)

new_field.SetWidth(32)
new_field.SetPrecision(2) #added line to set precision
layer.CreateField(new_field)

for feature in layer:
    geom = feature.GetGeometryRef()
    area = geom.GetArea() 
    print area
    feature.SetField("Area", area)
    layer.SetFeature(feature)


dataSource = None 

and it worked (see next image).

However, the precision of values (0 decimal) at the field "Area" is different to values printed at the Python Console:

1062218109.64
1241319130.43 

As you are pointed out that your printed areas are very small (0.00000x) and likely do not reflect square meters, this is the reason for your resulting "Area" field contains all 0's. Probably, you have a projection problem in your shapefile. It is not in meters.

Editing Note:

I included the code line to set the precision (2 decimals) of 'Area' field and it worked.