PyQGIS and the Python Data Stack

QGIS has its own complete Python API, and for a long time scripts that used it lived in a world apart from the rest of Python's data tools. That separation no longer makes sense. GeoPandas, shapely, pandas, NumPy, matplotlib and rasterio are the everyday tools of spatial data science, and they install into the same Python that runs QGIS. A script can load and clean data with PyQGIS, hand it to pandas for a pivot table, to shapely for vectorised geometry, to scikit-learn for a model, and bring the results back to QGIS to style and print — all in one process.

This guide belongs to Spatial Data Processing & Automation. It is for QGIS users who know some pandas or GeoPandas, and for Python data people who have started using QGIS. It explains how data crosses between the two worlds without losing CRS, NULLs, types or precision, and when each side is the better place for a given step.

QGIS objects and their Python counterpartsEach QGIS object has a natural counterpart in the Python data stack. A vector layer corresponds to a GeoDataFrame; its attribute table without geometry to a pandas DataFrame; a QgsGeometry to a shapely geometry; a raster layer's band to a NumPy array read with rasterio; and QGIS charts and layouts to matplotlib figures. WKB, file paths and feature ids are the bridges that carry data across.Each QGIS object has a counterpartQgsVectorLayerattribute tableQgsGeometryQgsRasterLayer bandGeoDataFramepandas DataFrameshapely geometryNumPy array (rasterio)path or WKBfeature id indexWKBfile path + profileresults come back the same way — by file, by WKB, by feature id

What this guide covers

One Python, two worlds

Everything starts with getting the libraries into the right interpreter. QGIS ships its own Python; packages installed into a system Python or a separate virtual environment are invisible to it. How to install depends on how QGIS was installed: the OSGeo4W shell on Windows, pip against QGIS's interpreter on macOS, the system package manager or pip on Linux, and conda install into the QGIS environment on conda-based setups. Installing Python packages into QGIS has the details per platform.

Two packages need extra care because they contain compiled GDAL or GEOS code: rasterio and, to a lesser extent, GeoPandas's I/O engines. They must be built against a GDAL compatible with the one QGIS uses, or importing them can crash QGIS. On conda and Linux distribution installs this is automatic when everything comes from the same channel; on Windows, use the OSGeo4W packages where they exist.

import sys
for name in ("numpy", "pandas", "shapely", "geopandas", "matplotlib", "rasterio"):
    try:
        mod = __import__(name)
        print(f"{name:<11} {mod.__version__}")
    except ImportError:
        print(f"{name:<11} not installed")
print("python", sys.version.split()[0], "at", sys.executable)

Breakdown: Running this in the QGIS Python console shows exactly which libraries the QGIS interpreter can see and where that interpreter lives — the path to use with pip if you install that way. NumPy is always present because QGIS depends on it; matplotlib and pandas are included in most official installers; the rest depend on the platform. For scripts outside QGIS Desktop that use the same libraries, see running Python scripts outside QGIS Desktop.

The three bridges

Data crosses between QGIS and the data stack over three bridges, and almost every conversion uses one of them.

Files. A layer backed by a GeoPackage, shapefile or GeoTIFF can be read by GeoPandas or rasterio directly from its path, and results written by them can be loaded by QGIS. This is the fastest route for large data and the most robust, because both sides use GDAL. The path comes from the layer's source string, decoded with QgsProviderRegistry.decodeUri. Its limitation is that it sees only what is on disk — not unsaved edits, subset filters, joins or virtual fields.

WKB. Geometry crosses between QgsGeometry and shapely as well-known binary, which both read and write natively and losslessly. Collecting WKB for a whole layer and converting it in one vectorised call is the efficient pattern.

Feature ids. Attributes cross as rows, and results come back by feature id. Indexing a DataFrame by the id each row came from is the single most important design decision in mixed scripts: it makes writing computed columns back a one-line dictionary comprehension, and it survives filtering, sorting and merging on the pandas side.

Choosing a bridgeFiles: fastest for large data, both sides use GDAL, but only sees saved data without subsets or joins. WKB: lossless geometry, vectorised, works for memory layers and edits, but attributes travel separately. Feature ids: carry results back to exactly the right features, survive filtering and merging, and allow one-call writes.Files, WKB, feature idsfilesfastest, GDAL both sidessees disk only:no edits, joins, subsetsWKBlossless geometryvectorisedcurves segmentized firstfeature idsresults back tothe right featuresone-call writes

What can go wrong at the border

Conversions fail quietly in a few predictable ways. Knowing them is most of the skill.

The CRS gets lost. WKB never carries a CRS, and a GeoDataFrame built without one has crs=None. Every later reprojection or overlay then misbehaves. Carry the layer's CRS across explicitly — as an EPSG code where it has one, as WKT otherwise — and check gdf.crs after every conversion.

NULLs and dates arrive as Qt objects. On QGIS 3, empty attributes are NULL QVariants and dates are QDate/QDateTime. Pandas stores them in object-dtype columns where numeric operations fail or, worse, silently skip. Clean them once at the border into None and Python datetimes, then let convert_dtypes() choose nullable pandas types. The rules are set out in handling NULL values and QVariant.

Types widen or narrow on the way back. An integer column with missing values becomes float in pandas; a 64-bit id written to a 32-bit field overflows; a timestamp becomes text. Map dtypes to field types deliberately when building layers, and verify counts, extents and a sample of values after loading.

Geometry types mix. A GeoDataFrame can hold points and polygons in one column; a layer cannot. Overlays in GeoPandas often produce stray lines and points where shapes touch. Split by geometry family, or filter, before loading.

Curves and M values disappear. Shapely has no circular arcs and limited M support. Segmentize curves in QGIS before converting, with a tolerance tied to the data's accuracy.

Which side should do the work?

With both toolkits available, the question for each step is where it runs best. A few rules of thumb keep scripts fast and simple.

Dividing the workQGIS side: loading from providers including databases and web services, standard geoprocessing on large layers through Processing, rendering, styling, layouts and atlases, and anything users will reproduce in the GUI. Python data stack side: group-bys, pivots and merges on attributes, vectorised measures over many geometries, statistical models and machine learning, array arithmetic on rasters, and charts. Convert once at the boundary between the two.Convert once, at the boundaryQGISload from any providerProcessing on big layersstyling, layouts, atlassteps users repeat in the GUIdata stackgroup-by, pivot, mergevectorised geometrymodels, arrays, chartssteps that are pure analysis

Keep loading and standard geoprocessing in QGIS. Providers read databases, web services and dozens of formats with authentication and caching handled; Processing algorithms run buffers, dissolves and overlays in compiled code and record parameters in the history. Converting to GeoPandas just to buffer gains nothing.

Move tabular and statistical work to pandas. Group-bys, pivots, merges with spreadsheets, rolling windows and period comparisons are one-liners in pandas and laborious anywhere else.

Use shapely when vectorisation pays. Measures, predicates and bulk spatial queries over tens of thousands of geometries are much faster as array calls than as Python loops over QgsGeometry.

Use NumPy and rasterio for raster arithmetic that does not fit an expression. Conditional masks, lookup-table reclassification, model predictions per pixel and block-wise processing of very large rasters.

Bring results back to QGIS for anything visual. Styling, labelling, layouts and atlas exports are QGIS's strengths, and matplotlib charts drop straight into layouts as SVG.

The pattern that follows is to convert once at each boundary rather than shuttle data back and forth between every step. A typical analysis has one conversion into the data stack and one back.

A complete round trip

The pieces fit together into short scripts. This one reads a layer, computes a per-district summary and a per-feature score with pandas, writes the score back, and saves a chart for the report.

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import pandas as pd
from qgis.core import QgsProject, QgsFeatureRequest, QgsField
from qgis.PyQt.QtCore import QVariant

layer = QgsProject.instance().mapLayersByName("buildings")[0]
names = ["district", "floor_area_m2", "energy_kwh"]
req = (QgsFeatureRequest().setFlags(QgsFeatureRequest.NoGeometry)
       .setSubsetOfAttributes(names, layer.fields()))
rows = {f.id(): [f[n] for n in names] for f in layer.getFeatures(req)}
df = pd.DataFrame.from_dict(rows, orient="index", columns=names).apply(
    pd.to_numeric, errors="ignore")

df["intensity"] = df["energy_kwh"] / df["floor_area_m2"]
district_median = df.groupby("district")["intensity"].transform("median")
df["vs_district"] = df["intensity"] / district_median

prov = layer.dataProvider()
if layer.fields().indexOf("vs_district") < 0:
    prov.addAttributes([QgsField("vs_district", QVariant.Double)])
    layer.updateFields()
i = layer.fields().indexOf("vs_district")
prov.changeAttributeValues({int(fid): {i: None if pd.isna(v) else float(v)}
                            for fid, v in df["vs_district"].items()})

ax = df.groupby("district")["intensity"].median().sort_values().plot.barh(color="#0f766e")
ax.set_xlabel("median kWh per m²")
plt.tight_layout()
plt.savefig("/data/reports/intensity_by_district.svg")
plt.close()

Breakdown: The dictionary keyed by feature id becomes the DataFrame index directly. transform("median") broadcasts each district's median back to its rows, so every building gets a ratio to its own district — the kind of calculation that is awkward in expressions and trivial in pandas. One provider call writes the ratio to every feature, ready for a graduated style that highlights buildings using far more energy than their neighbours. The chart is saved as SVG for a layout. Each step is covered in more depth in the recipes above.

Models: scikit-learn on layer data

Machine learning on spatial data follows the same bridges. Features — attributes plus geometric measures such as area, compactness or distance to the nearest road — go out as a NumPy matrix ordered by feature id; predictions come back as a column. QGIS supplies what the model cannot: spatial joins to build the predictors, and a map to see where the model is wrong.

import numpy as np
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score

parcels = QgsProject.instance().mapLayersByName("parcels_with_predictors")[0]
cols = ["area_m2", "dist_road_m", "slope_deg", "soil_index"]
data = {f.id(): ([f[c] for c in cols], f["land_use"]) for f in parcels.getFeatures()}
fids = [fid for fid, (x, y) in data.items() if y and all(v is not None for v in x)]
X = np.array([data[fid][0] for fid in fids], dtype=float)
y = np.array([data[fid][1] for fid in fids])

model = RandomForestClassifier(n_estimators=200, random_state=42)
print("cv accuracy", cross_val_score(model, X, y, cv=5).mean().round(3))
model.fit(X, y)

Breakdown: Building the matrix from a dictionary keyed by feature id keeps predictors and labels aligned, and the fids list records which rows made it through the NULL filter — the list you need to write predictions back. A fixed random_state makes the model reproducible. Ordinary cross-validation overstates accuracy on spatial data, because neighbouring parcels are similar and leak information between folds; group folds by district or a spatial grid for an honest estimate. Mapping the residuals — where predictions are wrong — usually teaches more than the accuracy figure.

Databases: SQL straight into pandas

When data lives in PostGIS, the shortest path to a DataFrame skips QGIS layers altogether: run SQL with the database driver and read the result directly. QGIS still helps by supplying the connection details from a stored connection, so credentials are not written into scripts.

from qgis.core import QgsProviderRegistry

md = QgsProviderRegistry.instance().providerMetadata("postgres")
conn = md.findConnection("city_db")
uri = conn.uri()
rows = conn.executeSql("SELECT district, count(*) AS n FROM permits "
                       "WHERE issued >= '2026-01-01' GROUP BY district")
permits = pd.DataFrame(rows, columns=["district", "n"])
print(permits.sort_values("n", ascending=False).head())

Breakdown: findConnection returns a connection saved in QGIS's browser, and executeSql runs a query and returns rows as lists — enough to build a DataFrame without installing a separate database driver. Aggregating in SQL means only the summary crosses the network. For geometry results, select ST_AsBinary(geom) and convert the bytes with shapely, the same WKB bridge as for layers. The connections API is covered in using the database connections API.

Notebooks for exploration

Exploratory work — trying a model, checking a distribution, iterating on a chart — is often more pleasant in a Jupyter notebook than in the QGIS console. PyQGIS runs in notebooks once the QGIS application is initialised, giving the same layers, providers and Processing algorithms alongside inline charts and tables. Using PyQGIS in a Jupyter notebook shows the setup. A common rhythm is to explore in a notebook, then move the finished steps into a script or a Processing algorithm that runs unattended.

Performance notes

Mixed scripts are usually fast enough without thought, but large layers expose a few hot spots.

  • Iterating features is the slow part. Read only needed columns, skip geometry when not needed, and filter in the request so the provider does the work.
  • Convert in batches. One shapely.from_wkb over a list beats a per-feature conversion by an order of magnitude; one changeAttributeValues beats a loop of single updates by more.
  • Prefer the file route for big file-backed layers. geopandas.read_file and rasterio block reads run in compiled code end to end.
  • Mind memory. A GeoDataFrame of a few million polygons can take gigabytes. Filter at the source, or process in chunks with keyset paging, as in reading feature attributes and geometry.
  • Run long work off the main thread. In a plugin, wrap the conversion and analysis in a QgsTask and apply results to layers when it finishes.

Reproducibility

Scripts that mix libraries depend on more versions. Record them: print the versions of QGIS, Python and each library at the start of a run, or write them into the output's metadata. Pin versions in a requirements or environment file for scheduled jobs. Seed any randomness — sampling, k-means initialisation, train/test splits — so a rerun gives the same answer. These habits matter most precisely when results are going into reports other people will rely on.

import json, platform
from qgis.core import Qgis

def environment_record():
    record = {"qgis": Qgis.version(), "python": platform.python_version()}
    for name in ("numpy", "pandas", "shapely", "geopandas", "matplotlib", "rasterio", "sklearn"):
        try:
            record[name] = __import__(name).__version__
        except ImportError:
            pass
    return record

with open("/data/results/run_environment.json", "w") as fh:
    json.dump(environment_record(), fh, indent=2)

Breakdown: A small JSON file next to the outputs answers the question that always comes up months later — which versions produced this? Qgis.version() gives the full QGIS version string, and the loop records whichever libraries are installed without failing on the ones that are not. Writing the same record into a GeoPackage's metadata or a layout's footer works just as well.

Key takeaways

  • Install the data stack into QGIS's own Python, and keep GDAL-linked packages from the same distribution as QGIS.
  • Files, WKB and feature ids are the three bridges; choose the route per layer.
  • Carry the CRS explicitly, clean NULLs and Qt dates at the border, and map types deliberately on the way back.
  • Keep loading, geoprocessing and visual output in QGIS; move tables, statistics, vectorised geometry and arrays to the data stack.
  • Convert once at each boundary, in batches, and write results back by feature id in one call.
  • Record versions and seed randomness so mixed scripts stay reproducible.

Frequently Asked Questions

Can I use GeoPandas instead of PyQGIS entirely? For pure analysis on files, often yes. PyQGIS adds providers for databases and web services, Processing's algorithm library, styling, layouts and plugins — use it where those matter.

Will these libraries work in QGIS 4? Yes. They are independent of QGIS's Qt version; only the Qt-facing parts of a script, such as field types and embedded matplotlib widgets, need the QGIS 4 adjustments noted in each recipe.

Is there a ready-made converter between layers and GeoDataFrames? Not in QGIS core. The short helper functions in the recipes are the common approach and make every choice visible.

What about polars, DuckDB or xarray? They fit the same pattern: attributes out by feature id, geometry as WKB (DuckDB's spatial extension reads it), rasters as arrays. The bridges do not change.

Do plugins have to bundle these libraries? Plugins should declare them rather than bundle compiled packages; see bundling third-party dependencies for the options.