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]]
Best Answer
I ran your script (slightly modified) at the Python Console of QGIS:
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:
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.