PyQGIS Script to Calculate Polygon Areas

To calculate polygon areas with PyQGIS, use the QgsGeometry.area() method inside a feature iteration loop. The script below adds a numeric field to your layer, computes planar area in the layer's native coordinate reference system (CRS), and writes values directly to the attribute table. For accurate linear measurements, your dataset must use a projected CRS (e.g., UTM or a national grid). Geographic CRS layers (like WGS84/EPSG:4326) will return meaningless square degrees from geom.area(). This workflow is a foundational step in Vector Data Manipulation for automated GIS tasks.

Working Script

Run this in the QGIS Python Console (Plugins > Python Console) or save as a .py file and execute via the Processing Toolbox.

from qgis.core import QgsProject, QgsField, edit
from qgis.PyQt.QtCore import QVariant

# 1. Reference target polygon layer
layer_name = "your_polygon_layer"
layers = QgsProject.instance().mapLayersByName(layer_name)
if not layers:
    raise RuntimeError(f"Layer '{layer_name}' not found in project.")
layer = layers[0]

# 2. Add output field
field_name = "area_sqm"
if field_name not in layer.fields().names():
    layer.dataProvider().addAttributes([QgsField(field_name, QVariant.Double)])
    layer.updateFields()

# 3. Calculate and write areas
with edit(layer):
    for feat in layer.getFeatures():
        geom = feat.geometry()
        if geom.isGeosValid() and not geom.isEmpty():
            feat.setAttribute(field_name, geom.area())
            layer.updateFeature(feat)

print(f"Area calculation complete for {layer.featureCount()} features.")

How It Works

The script follows a standard reference → schema update → iteration → commit pattern. QgsGeometry.area() returns planar area based on the layer's CRS units (meters → square meters for projected CRS). The with edit(layer): context manager batches all changes into a single transaction, drastically improving performance over per-feature disk writes. The isGeosValid() check prevents runtime crashes from self-intersecting rings, duplicate vertices, or collapsed geometries. Multi-polygon features are automatically aggregated by GEOS, requiring no manual flattening.

Compatibility & Fallbacks

  • QGIS 3.x Required: Uses modern edit() context manager and QgsField API. Legacy QGIS 2.x startEditing()/commitChanges() patterns are deprecated.
  • CRS Dependency: geom.area() is strictly planar. For geographic CRS (EPSG:4326), results will be in square degrees. Reproject first, or use the geodesic fallback below.
  • Provider Support: Works with Shapefiles, GeoPackage, PostGIS, and FileGDB. Read-only sources will fail on addAttributes().

1. Geodesic Area Fallback (No Reprojection Needed) When you cannot change the layer CRS, switch to ellipsoidal calculations using QgsDistanceArea:

from qgis.core import QgsDistanceArea, QgsProject, edit

da = QgsDistanceArea()
da.setSourceCrs(layer.crs(), QgsProject.instance().transformContext())
da.setEllipsoid("WGS84")

with edit(layer):
    for feat in layer.getFeatures():
        if feat.geometry().isGeosValid():
            feat.setAttribute(field_name, da.measureArea(feat.geometry()))
            layer.updateFeature(feat)

2. Processing Algorithm Alternative If direct PyQGIS editing fails due to provider locks or permission issues, use the native field calculator algorithm. This approach is safer for batch workflows within Spatial Data Processing & Automation pipelines:

import processing
from qgis.core import QgsProject

result = processing.run("native:fieldcalculator", {
    'INPUT': layer,
    'FIELD_NAME': 'area_sqm',
    'FIELD_TYPE': 0,        # Float
    'FIELD_LENGTH': 20,
    'FIELD_PRECISION': 6,
    'FORMULA': '$area',
    'OUTPUT': 'TEMPORARY_OUTPUT',
})
QgsProject.instance().addMapLayer(result['OUTPUT'])

Troubleshooting

  • AttributeError: 'QgsVectorLayer' object has no attribute 'dataProvider': The layer is likely not fully loaded or is a raster. Verify with from qgis.core import QgsMapLayer; layer.type() == QgsMapLayer.VectorLayer.
  • QgsVectorLayerError: Could not add field: The data source is read-only or lacks write permissions. Export to GeoPackage first using QgsVectorFileWriter.
  • NaN or 0.0 results: Usually caused by invalid topologies. Run processing.run("native:fixgeometries", {'INPUT': layer, 'OUTPUT': 'memory:'}) before calculating.
  • Always verify results by sampling polygons with known dimensions or comparing against the QGIS Identify tool. Wrap production runs in try/except and log failed feature IDs to maintain data integrity.

Frequently Asked Questions

Why does geom.area() return tiny decimal numbers for my WGS84 layer?QgsGeometry.area() is strictly planar and returns area in the layer's native CRS units. For a geographic CRS like EPSG:4326 those units are degrees, so you get square degrees—a meaningless value for real-world area. Reproject the layer to a projected CRS such as UTM first, or use the QgsDistanceArea geodesic fallback shown above.

What units does the calculated area come out in? With the planar geom.area() method the units follow the CRS: a metric projected CRS yields square meters, while a CRS in feet yields square feet. To convert to hectares, divide square meters by 10,000; for square kilometers, divide by 1,000,000. The QgsDistanceArea approach returns square meters when you set a metric ellipsoid like WGS84.

How do I get accurate areas without changing my layer's CRS? Use QgsDistanceArea with setSourceCrs() and setEllipsoid("WGS84"), then call da.measureArea(geom). This performs an ellipsoidal (geodesic) calculation that stays accurate regardless of the layer's stored projection, so you never have to physically reproject the data.

Why are some of my area values 0.0 or NaN? These results almost always come from invalid topologies—self-intersecting rings, collapsed geometries, or empty features. Run processing.run("native:fixgeometries", {'INPUT': layer, 'OUTPUT': 'memory:'}) before calculating, and keep the isGeosValid() and isEmpty() guards in the loop so problem features are skipped rather than corrupting the field.

Can I calculate areas without writing a script? Yes. Use the native field calculator algorithm with the $area expression, either through processing.run("native:fieldcalculator", ...) as shown above or via the Processing Toolbox GUI. This is often safer for batch workflows because it avoids provider locks and writes to a fresh output layer.