Short answer

Python GIS and desktop GIS solve the same spatial problems through opposite interfaces. Desktop GIS (QGIS, ArcGIS Pro) is interactive and visual: you see every layer and every intermediate result on a canvas. Python GIS (GeoPandas, Shapely, Rasterio, PyProj) is code: no map, but full reproducibility, batch scale, and integration into pipelines. You script when an analysis must be rerun, applied across many files, version-controlled, or embedded in a larger system; you stay on the desktop for exploration, editing, cartography, and quick one-offs. The professional norm is to use both, and the two are not even mutually exclusive — both desktop tools run Python internally.

Two interfaces to the same concepts

Nothing conceptual differs between the two. A coordinate reference system, a geometry, a spatial join, a buffer, a raster band — these are identical whether you click them or call them. What differs is the loop.

In desktop GIS the loop is see, click, see. You load a layer, watch it draw, run a tool, inspect the output visually, adjust. This is unbeatable when you are still figuring out what the data is or what you want, and it is the right environment for anything that needs human judgement on each result: digitising, fixing topology by hand, designing a print layout, tuning symbology.

In Python the loop is write, run, verify. You describe the whole operation in code, execute it, and check the result programmatically. There is no canvas to glance at (you plot deliberately), which makes the start slower — but the code is a permanent, exact, shareable record of what was done.

What the Python stack gives you

The core libraries map cleanly onto desktop functionality:

  • GeoPandas — vector data as a table with a geometry column. Reads/writes shapefiles, GeoPackages, GeoJSON, PostGIS; does spatial joins, dissolves, overlays, reprojection. It is "the attribute table plus the processing toolbox" as a DataFrame.
  • Shapely — geometry operations on individual shapes: buffer, intersection, union, distance, validity.
  • Rasterio — raster I/O and array access, built on GDAL; read windows, reproject, compute on bands as NumPy arrays.
  • PyProj — CRS definitions and coordinate transformations (the PROJ library).
  • GDAL/OGR and Fiona/pyogrio — format conversion and low-level access.

A vector workflow in GeoPandas reads almost like the QGIS steps it replaces:

import geopandas as gpd

parcels = gpd.read_file("parcels.gpkg")
faults  = gpd.read_file("faults.gpkg")

# put both in a metric CRS so distances are in metres
parcels = parcels.to_crs(25832)
faults  = faults.to_crs(25832)

# parcels within 500 m of a fault
near = gpd.sjoin_nearest(parcels, faults, max_distance=500, how="inner")
near.to_file("hazard_parcels.gpkg")

That is the same buffer-and-select-by-location task a QGIS user does by hand — but it is now a file you can rerun, diff, review, and schedule. Note the same CRS discipline applies in both worlds: to_crs(25832) is essential because max_distance=500 is only metres in a projected CRS.

When to script, when to click

Reach for Python when:

  • The same operation must run repeatedly (monthly updates, audited deliverables).
  • You must process many files in batch — reprojecting 400 DEM tiles is a for loop, not 400 click-throughs.
  • The result must be reproducible and reviewable as code.
  • The GIS step is one stage of a bigger pipeline (ingest → process → publish), or must run on a server with no display.
  • The dataset is large and you want NumPy/array-level control over rasters.

Stay in desktop GIS when:

  • You are exploring unfamiliar data and need to see it.
  • The task needs interactive editing or hand digitising.
  • You are making a map for humans — layout, labels, symbology.
  • It is a genuine one-off where writing code costs more than it saves.

The honest test is recurrence and judgement: recurring, mechanical work belongs in code; one-off, judgement-heavy work belongs on the canvas.

They are not separate worlds

You do not have to choose at the boundary. QGIS embeds a full Python console and the PyQGIS API, and every Processing algorithm is callable from Python (processing.run("native:buffer", {...})). ArcGIS Pro exposes the same via ArcPy. So the natural workflow is to prototype interactively on the canvas, then lift the working steps into a script for repeatability — you get the visual feedback of the desktop and the automation of code from the same toolset. Many teams develop a Processing model visually, then export it to a Python script.

A worked batch example

Reproject and compute slope for every DEM tile in a directory — a task that is tedious in the GUI and trivial in code:

import glob, subprocess, os

for src in glob.glob("dem_tiles/*.tif"):
    base = os.path.splitext(os.path.basename(src))[0]
    proj = f"proj/{base}.tif"
    subprocess.run(["gdalwarp", "-t_srs", "EPSG:25832",
                    "-tr", "10", "10", "-r", "bilinear", src, proj], check=True)
    subprocess.run(["gdaldem", "slope", proj,
                    f"slope/{base}.tif", "-compute_edges"], check=True)

This wraps the same GDAL commands a desktop user would run via the GUI, but applies them uniformly and reproducibly across the whole tile set. Scale is exactly where scripting earns its keep.

Common pitfalls and why they happen

  • Forcing everything into code. Scripting an exploratory, one-time look at unfamiliar data is slower than just opening QGIS; the value of code is reuse, and a one-off has none.
  • Refusing to script recurring work. A monthly process done by hand drifts — slightly different parameters each time — producing inconsistent results. That is a reproducibility failure disguised as familiarity.
  • Ignoring CRS in code. Without to_crs(), GeoPandas/Shapely happily compute distances in degrees. The library will not warn you; the canvas at least shows misalignment.
  • Reinventing what GDAL already does. Many "Python GIS" tasks are best as a subprocess call to a battle-tested GDAL program rather than hand-rolled raster math.

Validation

  • Run the same operation in QGIS and in Python on a sample and confirm matching feature counts, extents, and a few computed values — a fast catch for CRS or geometry-validity differences.
  • Assert CRS explicitly in scripts (assert gdf.crs.to_epsg() == 25832) rather than trusting inputs.
  • Keep analysis scripts in version control; the script is the audit trail.
  • Plot or write a sample output and eyeball it — scripting removes the canvas, so build a deliberate visual check back in.

Bathyl perspective

We script when GIS work becomes a system — recurring, batched, or part of a pipeline — and keep the desktop for the parts that genuinely need a human looking at a map. Because QGIS and ArcGIS both run Python natively, the boundary is soft: prototype on the canvas, automate in code, and let the script be the durable record so the same result can be reproduced by anyone, on next month's data.

Related reading

Sources