Automating QGIS with PyQGIS

"""A standalone PyQGIS script: QGIS started without a window, a GeoPackage read, and QGIS closed again."""
import os
import sys

os.environ.setdefault("QT_QPA_PLATFORM", "offscreen")  # draw without a screen, as on a server
os.environ.setdefault("QT_LOGGING_RULES", "qt.*=false")  # keep Qt's own notes out of the output
try:
    from qgis.core import QgsApplication
except ImportError:  # a conda-forge QGIS keeps its Python package in share/qgis/python
    sys.path.insert(0, os.path.join(sys.prefix, "share", "qgis", "python"))
    from qgis.core import QgsApplication

QgsApplication.setPrefixPath(os.environ.get("QGIS_PREFIX_PATH", sys.prefix), True)  # where QGIS is installed
app = QgsApplication([], False)  # no command-line arguments, and False: no window
app.initQgis()  # loads the data providers and the layer registry

from qgis.core import Qgis, QgsProviderRegistry, QgsVectorLayer, QgsWkbTypes  # noqa: E402

from town_files import write_town  # noqa: E402

print("QGIS", Qgis.version().split("-")[0], "| data providers:", len(QgsProviderRegistry.instance().providerList()))
write_town("town.gpkg")
for name in ("wards", "households", "clinics"):
    layer = QgsVectorLayer(f"town.gpkg|layername={name}", name, "ogr")
    print(f"  {name:10} valid {layer.isValid()}, {layer.featureCount():5,} features,",
          f"{QgsWkbTypes.displayString(layer.wkbType())}, {layer.crs().authid()}")
    del layer  # the layer goes before QGIS does
app.exitQgis()  # lets go of the data providers before Python exits
print("exitQgis: done")

QGIS 4.2.3 | data providers: 34
  wards      valid True,    12 features, Polygon, EPSG:32643
  households valid True, 5,000 features, Point, EPSG:32643
  clinics    valid True,     8 features, Point, EPSG:32643
exitQgis: done

"""Five feature requests: ask a layer for only the features, the fields and the geometry that a task needs."""
import qgis_start
from qgis.core import Qgis, QgsExpression, QgsFeatureRequest, QgsGeometry, QgsVectorLayer

from town_files import write_town

write_town("town.gpkg")
households = QgsVectorLayer("town.gpkg|layername=households", "households", "ogr")
wards = QgsVectorLayer("town.gpkg|layername=wards", "wards", "ogr")

# 1. An attribute filter: an expression, evaluated by the data provider where it can be.
w04 = next(wards.getFeatures(QgsFeatureRequest(QgsExpression("\"ward\" = 'W04'"))))
print("Ward W04 found by an expression; its box:", w04.geometry().boundingBox().toString(0))

# 2. A spatial filter: features whose box meets a rectangle, found through the GeoPackage's spatial index.
request = QgsFeatureRequest().setFilterRect(w04.geometry().boundingBox())
inside = sum(1 for _ in households.getFeatures(request))
everything = sum(1 for _ in households.getFeatures())
print(f"Households in W04's box: {inside:,}, read instead of all {everything:,}")

# 3. A rectangle stays a rectangle: an L of three wards, W03, W04 and W08, has a fourth ward, W07, inside its box.
shapes = {f["ward"]: f.geometry() for f in wards.getFeatures()}
l_shape = QgsGeometry.unaryUnion([shapes["W03"], shapes["W04"], shapes["W08"]])
box_only = QgsFeatureRequest().setFilterRect(l_shape.boundingBox())
exact = QgsFeatureRequest(box_only).setFlags(Qgis.FeatureRequestFlag.ExactIntersect)  # exact, but against the box
in_l = QgsFeatureRequest().setDistanceWithin(l_shape, 0)  # the index finds candidates, the L's own shape decides
counts = [sum(1 for _ in households.getFeatures(request)) for request in (box_only, exact, in_l)]
print("An L of W03, W04 and W08: {:,} households in its box, {:,} with ExactIntersect, {:,} inside the L".format(*counts))

# 4. Only some fields, and no geometry at all, when the task reads one attribute.
names = QgsFeatureRequest()
names.setSubsetOfAttributes(["household"], households.fields())
names.setFlags(Qgis.FeatureRequestFlag.NoGeometry)
listed = [f for f in households.getFeatures(names)]
print(f"NoGeometry and one field: {len(listed):,} names read, first {listed[0]['household']},",
      "any geometry read:", any(f.hasGeometry() for f in listed))
both = QgsFeatureRequest(names).setFilterRect(w04.geometry().boundingBox())
in_box = [f for f in households.getFeatures(both)]
print(f"The same request inside W04's box: {len(in_box):,} names; geometry read for the box test:",
      all(f.hasGeometry() for f in in_box))

# 5. A limit: the first few features, for a preview or a quick test.
preview = [f["household"] for f in households.getFeatures(QgsFeatureRequest().setLimit(3))]
print("setLimit(3):", ", ".join(preview))

Ward W04 found by an expression; its box: 378000,2048000 : 380000,2050000
Households in W04's box: 1,169, read instead of all 5,000
An L of W03, W04 and W08: 2,071 households in its box, 2,071 with ExactIntersect, 1,692 inside the L
NoGeometry and one field: 5,000 names read, first H0001, any geometry read: False
The same request inside W04's box: 1,169 names; geometry read for the box test: True
setLimit(3): H0001, H0002, H0003

"""Editing a layer in a session: changes that are kept together on success, or all thrown away on an error."""
import qgis_start
from qgis.core import QgsField, QgsVectorLayer, edit
from qgis.PyQt.QtCore import QMetaType

from town_files import write_town

write_town("town.gpkg")
wards = QgsVectorLayer("town.gpkg|layername=wards", "wards", "ogr")


def saved_areas():
    """The area_km2 values as the file holds them, read by a fresh layer."""
    fresh = QgsVectorLayer("town.gpkg|layername=wards", "check", "ogr")
    return [f["area_km2"] for f in fresh.getFeatures()][:3]


with edit(wards):  # starts an edit session, and commits it when the block ends without an error
    wards.addAttribute(QgsField("area_km2", QMetaType.Type.Double))
    wards.updateFields()
    column = wards.fields().indexOf("area_km2")
    for ward in wards.getFeatures():
        wards.changeAttributeValue(ward.id(), column, round(ward.geometry().area() / 1e6, 2))
print("Committed: the file now holds area_km2 =", saved_areas(), "...")

try:
    with edit(wards):
        for ward in wards.getFeatures():
            wards.changeAttributeValue(ward.id(), column, 0.0)
            if ward["ward"] == "W02":
                raise ValueError("a check failed at W02, halfway through")
except ValueError as error:
    print("The second session stopped:", error)
print("Rolled back: still editing?", wards.isEditable(), "| the file holds", saved_areas(), "...")

# The data provider writes at once: no session, no undo.
wards.dataProvider().changeAttributeValues({next(wards.getFeatures()).id(): {column: 9.99}})
print("Written straight through the data provider:", saved_areas(), "...")

Committed: the file now holds area_km2 = [4.0, 4.0, 4.0] ...
The second session stopped: a check failed at W02, halfway through
Rolled back: still editing? False | the file holds [4.0, 4.0, 4.0] ...
Written straight through the data provider: [9.99, 4.0, 4.0] ...

"""The toolbox's algorithms called from Python: the same ids and parameter names, and a dictionary of results back."""
import qgis_start
from qgis.core import QgsVectorLayer

from town_files import write_town

processing = qgis_start.processing()
write_town("town.gpkg")

near = processing.run("native:extractwithindistance", {
    "INPUT": "town.gpkg|layername=households", "REFERENCE": "town.gpkg|layername=clinics", "DISTANCE": 1000,
    "OUTPUT": "memory:"})["OUTPUT"]
print("Extract within distance returned a layer of", f"{near.featureCount():,}", "households")

result = processing.run("native:countpointsinpolygon", {
    "POLYGONS": "town.gpkg|layername=wards", "POINTS": near, "FIELD": "near", "OUTPUT": "near.gpkg"})
print("Count points in polygon returned", result)
counted = QgsVectorLayer(result["OUTPUT"], "near", "ogr")
print("Per ward:", ", ".join(f"{f['ward']} {f['near']:.0f}" for f in counted.getFeatures()))

Extract within distance returned a layer of 2,973 households
Count points in polygon returned {'OUTPUT': 'near.gpkg'}
Per ward: W01 210, W02 304, W03 100, W04 691, W05 27, W06 320, W07 313, W08 146, W09 94, W10 0, W11 444, W12 324