5

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]]

  • Do you get a result when you print geom.GetArea() in the for loop? – Master Bee Nov 5 '15 at 18:21
  • Yes! Areas print to the command line – geoeye Nov 5 '15 at 18:24
  • However the numbers are very small (0.00000x) and likely do not reflect square meters. – geoeye Nov 5 '15 at 18:41
  • 1
    I cannot test this right now, but try either setting a precision or not setting a width. Since your projection should not be the problem, maybe you are truncating your real with a very large precision? You can also try buffering the geometry by zero (in case of poor/invalid geometry) and printing the spatial ref of each geometry to make sure it is still your original shapefile projection. – blord-castillo Nov 5 '15 at 22:13
  • Removing SetWidth does the trick! I now get an area field with very small values. However, my projected coordinate system is definitely in units meters as demonstrated by info in the .prj file (see question). Any ideas why these values do not reflect meters? – geoeye Nov 6 '15 at 20:06
6

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

enter image description here

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.

  • Setting precision instead of width did the trick! I now get an area field appended on with very small values. However, my projected coordinate system is definitely in units meters as demonstrated by info in the .prj file (see question). Any idea why this discrepancy may arise? I tried reprojecting the shapefile before calculating the area field but to no avail. – geoeye Nov 6 '15 at 20:06
  • Info in the .prj file is not demonstrating that your coordinate system is definitely in meters. It is indicating an assigned projection (It could be erroneous). The facts (areas values) are indicating that the assigned shapefile projection seems not correct. – xunilk Nov 6 '15 at 21:00
  • 1
    How best to go about fixing this? – geoeye Nov 6 '15 at 21:08
  • What is the origin of your shapefile? Was it reprojected by you? – xunilk Nov 7 '15 at 0:08
  • Yes. It is a polygonized raster representative of flood extent. I reprojected it when polygonizing it using the t_srs argument in gdal.polygonize. – geoeye Nov 8 '15 at 20:56
1

Your code is perfect. I think you are not working in a projection in meters e.g mollweide.

All you need to do is reproject the layer to -- +proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m no_defs

Then run the code again and it would be 100%

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.