PyQGIS Script to Calculate Polygon Areas

To calculate polygon areas with PyQGIS, iterate the features of a polygon layer and read each geometry's area, then write the value into a numeric attribute field. The catch that trips up almost everyone is the coordinate reference system: QgsGeometry.area() returns planar area in the layer's own CRS units, so a projected CRS (UTM, a national grid) gives square meters while a geographic CRS such as WGS84/EPSG:4326 gives meaningless square degrees. This page shows a correct, runnable script, the geodesic alternative that stays accurate without reprojecting, and the pitfalls that produce 0.0 or NaN values.

Adding a computed area column is one of the most common tasks in Vector Data Manipulation: sizing parcels, reporting land-cover totals, or filtering polygons above a threshold. The same pattern — reference a layer, add a field, iterate, commit — underpins most attribute-writing scripts, so it is worth getting right once.

Prerequisites

  • QGIS 3.34 LTR (bundled Python 3.12), or any current 3.x release.
  • A polygon layer loaded in the project (Shapefile, GeoPackage, PostGIS, or FileGDB) with a writable data source.
  • Familiarity with the QGIS Python Console (Plugins > Python Console, or Ctrl+Alt+P).
  • Ideally, a projected CRS on the layer so planar area comes out in linear units. Check with layer.crs().isGeographic() — if it returns True, use the geodesic recipe further down instead of reprojecting.

In the Python Console iface and processing are already imported. In a standalone script you must first initialize QgsApplication and Processing — see Run a Processing Algorithm from a Script in PyQGIS.

Add an Area Field and Populate It

Run this in the Python Console, or save it as a .py file and execute it as a script. It references the target layer by name, adds a Double field if it is missing, then writes the planar area of every valid feature inside a single edit transaction.

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

# 1. Reference the 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 the output field if it does not already exist
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 inside one transaction
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.")

Breakdown: mapLayersByName resolves the layer by its tree name; guarding against an empty list gives a clear error instead of an IndexError. QgsGeometry.area() returns planar area in the CRS's units — square meters for a metric projected CRS. The with edit(layer): context manager batches every change into one commit, which is dramatically faster than per-feature disk writes and rolls back cleanly if an exception is raised. The isGeosValid() and isEmpty() guards skip self-intersecting rings, collapsed geometries, and empty features that would otherwise poison the result. Multipolygons are aggregated by GEOS automatically, so no manual flattening is needed.

Geodesic Area Without Reprojecting

When the layer is in a geographic CRS and you cannot (or do not want to) reproject it, switch from planar to ellipsoidal measurement with QgsDistanceArea. This is the accurate way to measure large polygons that span several UTM zones, and it is closely tied to correct coordinate reference system handling.

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():
        geom = feat.geometry()
        if geom.isGeosValid() and not geom.isEmpty():
            feat.setAttribute(field_name, da.measureArea(geom))
            layer.updateFeature(feat)

Breakdown: setSourceCrs tells the calculator what projection the incoming geometries use, and setEllipsoid("WGS84") selects the reference ellipsoid for the on-the-ellipsoid math. da.measureArea() then returns square meters regardless of the layer's stored CRS, so you never have to physically reproject the data. Newer QGIS builds may emit a deprecation note for measureArea; if so, da.measureAreaGeometry(geom) is the drop-in successor.

Planar versus geodesic area decision in PyQGISStarting from a polygon layer, test layer.crs().isGeographic(). If false the CRS is projected, so QgsGeometry.area() returns valid square metres. If true the CRS is geographic such as EPSG:4326, where planar area() returns meaningless square degrees; instead either reproject to UTM and call area(), or call QgsDistanceArea.measureArea() with the WGS84 ellipsoid, both of which return valid square metres.Polygon layerlayer.crs().isGeographic()?False · projectedTrue · geographicQgsGeometry.area()planar, in CRS linear units✓ Square metresfast & exact on a metric CRSPlanar area() here → square degrees✗ meaningless — degrees are not a lengthReproject to UTMthen call .area()✓ square metresQgsDistanceArea.measureArea() ·WGS84 ellipsoid✓ square metres

No-Script Alternative: the Field Calculator Algorithm

If direct editing fails because of a provider lock or a read-only source, the native field-calculator algorithm writes the area to a fresh output layer instead of mutating the original. This is often the safest option inside larger 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,        # 0 = Float
    "FIELD_LENGTH": 20,
    "FIELD_PRECISION": 6,
    "FORMULA": "$area",
    "OUTPUT": "TEMPORARY_OUTPUT",
})
QgsProject.instance().addMapLayer(result["OUTPUT"])

Breakdown: The $area expression uses the project's ellipsoid setting when one is configured, so it can already return ellipsoidal area without extra code. TEMPORARY_OUTPUT produces a scratch layer that leaves your source untouched — swap it for a .gpkg path to persist the result. This is the same approach used for other attribute jobs across the vector toolset, such as merging multiple shapefiles before summarizing them.

Convert to Hectares or Square Kilometers

Area rarely stays in square meters for reporting. Derive additional columns in the same loop rather than post-processing:

with edit(layer):
    for feat in layer.getFeatures():
        geom = feat.geometry()
        if geom.isGeosValid() and not geom.isEmpty():
            sqm = geom.area()
            feat.setAttribute("area_sqm", sqm)
            feat.setAttribute("area_ha", sqm / 10_000)      # hectares
            feat.setAttribute("area_sqkm", sqm / 1_000_000)  # square km
            layer.updateFeature(feat)

Breakdown: One square meter is 1e-4 hectares and 1e-6 square kilometers, so the conversions are plain division. Add the area_ha and area_sqkm fields the same way as area_sqm before the loop runs.

Planar area versus ellipsoidal area

geometry.area() measures on the flat plane of the layer's CRS. QgsDistanceArea with an ellipsoid set measures on the curved Earth. For a single parcel the two agree closely; for a region they do not, and for a layer in EPSG:4326 the planar figure is not an area at all.

Three area figures for one polygonOne parcel is measured three ways. In EPSG:4326 the area comes out as a tiny number in square degrees, which is not a unit of area. In a UTM zone the planar area is in square metres and close to correct. The ellipsoidal measurement is in square metres and accounts for the curvature of the Earth, differing from the planar figure by a fraction of a percent at parcel scale and much more at regional scale.Only two of these three numbers are areasEPSG:4326 · geometry.area()0.00000412square degrees — meaninglessEPSG:32633 · geometry.area()4 218.62 m²planar — good at parcel scaleQgsDistanceArea · measureArea()4 219.05 m²ellipsoidal — correct anywhere

The rule of thumb is simple: if the layer might ever be in a geographic CRS, or the polygons might ever be larger than a few square kilometres, use QgsDistanceArea. It works correctly on lat/lon data without reprojection, which alone usually justifies the three extra lines of setup.

Rounding, units and what to store

An area field is a derived value, and how it is stored decides whether it can be trusted a year later.

What to store alongside a computed areaStoring raw square metres at full precision is precise but implies accuracy the source geometry does not have. Storing rounded hectares alone loses the unit, so a later reader cannot tell what the number means. Storing a rounded value alongside an explicit units column and the calculation date is unambiguous and auditable.A number without its unit is a future support ticketraw m²4218.6249183precise to a micronthe survey was notbare number0.42hectares? acres?nobody can tell latervalue + unit + datearea_value 0.4219area_unit haarea_calc 2026-08-01unambiguous and auditable

Recording the calculation date matters more than it sounds. Areas go stale the moment geometry is edited, and a date column is the only cheap way to spot a value computed before the last round of boundary corrections.

QGIS Version Compatibility

The examples target QGIS 3.34 LTR (Python 3.12). The area APIs are stable across all current 3.x releases.

QGIS versionPythonNotes
3.28 LTR3.9edit(), QgsGeometry.area(), and QgsDistanceArea all identical.
3.34 LTR3.12Baseline for this page.
3.40 / 3.443.12Same API; QgsDistanceArea.measureArea() may warn — prefer measureAreaGeometry().

Two version notes matter in practice. First, the legacy QGIS 2.x startEditing() / commitChanges() pattern still works but is superseded by the edit() context manager shown here. Second, in QGIS 3.36+ QVariant.Double can be replaced by the newer QMetaType field constructor, but QVariant.Double remains valid and portable, so the code above runs unchanged on every LTR.

Troubleshooting

  • geom.area() returns tiny decimals (e.g. 0.00042). The layer is in a geographic CRS, so the result is square degrees. Reproject to a projected CRS, or use the QgsDistanceArea recipe above.
  • AttributeError: 'QgsVectorLayer' object has no attribute 'dataProvider'. The object is not a vector layer (often a raster or a partially loaded layer). Confirm with from qgis.core import QgsMapLayer; layer.type() == QgsMapLayer.VectorLayer.
  • Could not add field / write refused. The data source is read-only. Export to GeoPackage first with QgsVectorFileWriter, or use the field-calculator algorithm, which writes a new layer.
  • 0.0 or NaN values. Invalid topology — self-intersecting rings or collapsed geometries. Run processing.run("native:fixgeometries", {"INPUT": layer, "OUTPUT": "memory:"}) first, and keep the isGeosValid() / isEmpty() guards so bad features are skipped rather than corrupting the column.
  • Values look plausible but wrong. Verify against a polygon of known size or the QGIS Identify tool. In production, wrap the loop in try/except and log failing feature IDs so one bad geometry does not abort the whole run.

Conclusion

Calculating polygon areas in PyQGIS comes down to one decision: planar or geodesic. On a projected CRS, QgsGeometry.area() inside an edit() transaction is fast, exact, and only a few lines. On a geographic CRS, QgsDistanceArea with a WGS84 ellipsoid gives accurate square meters without touching the projection. Guard every geometry, batch writes into a single commit, and validate the output against a known reference — with those habits the same script scales cleanly from a single Console run to an automated batch pipeline.

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 such as EPSG:4326 those units are degrees, so you get square degrees — meaningless as real-world area. Reproject the layer to a projected CRS such as UTM, or use the QgsDistanceArea geodesic recipe shown above.

What units does the calculated area come out in? With planar geom.area() the units follow the CRS: a metric projected CRS yields square meters, a CRS in feet yields square feet. Divide square meters by 10,000 for hectares or by 1,000,000 for square kilometers. 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 ellipsoidal calculation 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 almost always come from invalid topology — self-intersecting rings, collapsed geometries, or empty features. Run native:fixgeometries before calculating and keep the isGeosValid() and isEmpty() guards in the loop so problem features are skipped rather than writing garbage into 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. It avoids provider locks by writing to a fresh output layer, which is often safer for batch workflows.