Spatial Processing on HPC: From Your Laptop to a Thousand Cores

You have been granted an account on a High Performance Computing (HPC) facility. Perhaps your institution bought into a national system, perhaps a collaborator has an allocation with spare hours. You log in, and every tutorial you find was written for someone running molecular dynamics.

This article is the geospatial version. It assumes you can write a Python script that processes a raster and that you have never used a job scheduler. By the end you will understand what a cluster actually is, why it will probably make your individual tasks slower, how the storage architecture changes the correct answer to nearly every design question, and how to take one working script from your laptop to five hundred cores.

# What an HPC Actually Is

An HPC facility is not a big computer. It is a few hundred or few thousand ordinary computers, wired together, with software deciding who gets to use which ones.

The pieces you need to know about:

Login nodes. Where your Secure Shell (SSH) session lands. Shared by every user on the system. You edit files here, submit jobs here, and you do not compute here. Running a heavy job on a login node is the fastest way to get an email from a system administrator.

Compute nodes. Where work happens. You never touch them directly; you ask the scheduler for some and it runs your script on them. Typically they have no outbound internet access, which surprises everyone the first time pip install hangs.

A parallel filesystem. Usually Lustre or the General Parallel File System (GPFS), presented as an ordinary directory like /scratch/ab1234. Every node sees the same files. It is built for streaming large files fast and, as we will see, has strong opinions about how you use it.

Node-local scratch. Some nodes have their own Non-Volatile Memory Express (NVMe) drive, exposed as $TMPDIR or similar. Extremely fast, private to your job, and erased the moment the job ends.

An interconnect. InfiniBand (IB), usually, rather than Ethernet. The relevant difference is latency — roughly one microsecond node-to-node versus fifty for Gigabit Ethernet (GbE). This is what makes tightly coupled multi-node work viable at all.

A module system. Lmod, a Lua-based implementation of Environment Modules, in most cases. Software is not installed system-wide; you load what you need with module load.

A scheduler. Slurm on most systems, sometimes Portable Batch System (PBS) Pro. It owns every compute node, maintains queues, enforces limits, and decides when your work runs. Everything you do on a cluster is mediated by it.

Time is charged against an allocation, usually denominated in core-hours or an abstract Service Unit (SU). Using 400 cores for an hour costs the same as 1 core for 400 hours.

# How It Differs From a Fast Desktop

This is where intuition misleads people, so here it is side by side.

Read the clock-speed row again, because it is the one that catches people. An HPC will usually make your job slower and your workflow harder. A single tile that takes 40 seconds on your desktop may take 80 on a compute node. The cluster wins only when the work decomposes into many independent pieces, because it can run five hundred of those pieces at once.

That gives a straightforward decision rule:

• Under roughly a thousand core-hours, and it fits on one machine — stay on a single fat box and use Dask-GeoPandas or multiprocessing. The queue wait alone will exceed your runtime.
• Embarrassingly parallel at large scale — thousands of tiles, scenes, or catchments with no interdependence. This is the sweet spot, and job arrays are all you need.
• Tightly coupled — the calculation for one region needs values from its neighbours at every step. This needs the Message Passing Interface (MPI) and is the case where the interconnect earns its cost.
• Bursty, needs the internet, or containerised — cloud batch is a better fit. See Cloud-Orchestrated Geospatial Workflows.

# Week One: Getting a Working Environment

Set up your connection first, because you will do it a hundred times:

# ~/.ssh/config
Host gadi
  HostName gadi.nci.org.au
  User ab1234
  IdentityFile ~/.ssh/id_ed25519
  ControlMaster auto
  ControlPath ~/.ssh/cm-%r@%h:%p
  ControlPersist 10m

Then orient yourself:

sinfo -s                       # partitions, and how busy they are
squeue -u $USER                # your jobs
sacctmgr show assoc user=$USER format=account,partition,grpsubmit
module avail gdal              # what geospatial software exists already
module load gdal/3.8.0

Now the geospatial-specific trap. A conda environment containing the Geospatial Data Abstraction Library (GDAL), GeoPandas, and Rasterio is roughly 150,000 small files. On a filesystem with a metadata bottleneck, importing from it across 400 concurrent tasks generates a metadata storm that degrades the filesystem for every other user on the machine. There are two acceptable answers.

Option A — a packed environment on fast storage. Build with micromamba, then collapse it:

micromamba create -y -p /scratch/ab1234/envs/geo -c conda-forge \
  python=3.11 gdal rasterio geopandas dask-jobqueue mpi4py
micromamba activate /scratch/ab1234/envs/geo
conda-pack -p /scratch/ab1234/envs/geo -o geo.tar.gz

Option B — a container. Docker is never permitted on a shared cluster, because the daemon runs as root. Apptainer (formerly Singularity) is the HPC-native equivalent and runs unprivileged, producing a single Singularity Image Format (SIF) file — one inode instead of 150,000:

apptainer build geo.sif docker://ghcr.io/osgeo/gdal:ubuntu-full-3.8.0
apptainer exec --bind /scratch:/scratch geo.sif python my_script.py

Prefer the container. It is reproducible, it is a single file, and it removes an entire class of “works for me” problems.

# Know Your Storage Tier Before You Write a Line of Code

Here is the part that most HPC guidance gets wrong by assuming one architecture. Storage is a latency ladder, and which rung your data sits on determines the correct design.

Note the shape of that table. The gap between tape and everything else is five or six orders of magnitude. Every other transition is one or two. The only step worth architecting around is the recall from tape. Once data is on flash, the chain of tape, staging, flash, memory, and processing collapses to “read it when you need it”.

# The Classic Case: Lustre or GPFS, No Flash Tier

On a traditional parallel filesystem, every open() and stat() is routed through a small number of Metadata Server (MDS) nodes, while the bulk data lives on many Object Storage Targets (OSTs). Bandwidth scales with the OSTs; metadata operations do not scale at all. Four hundred tasks each opening five thousand files will saturate the MDS regardless of how fast the underlying disks are.

So, on this architecture:

• Consolidate. One GDAL Virtual Format (VRT) file over 5,000 tiles instead of 5,000 separate opens. Cloud-Optimised GeoTIFF (COG), Zarr, and GeoParquet instead of directories of small files.
• Stripe large files across OSTs: lfs setstripe -c 8 /scratch/ab1234/bigdata/.
• Stage in, compute local, stage out. Copy inputs to $TMPDIR, work there, copy results back.
• Never point hundreds of ranks at one GeoPackage, SQLite database, or single PostGIS instance. File locking will serialise you or simply fail.
• One output file per rank, merged afterwards.
• Respect the memory ceiling. At 4 GB per core you use windowed reads, not src.read().
• Pin thread counts, or you will oversubscribe catastrophically: OMP_NUM_THREADS=1, GDAL_NUM_THREADS=1, and a sensible GDAL_CACHEMAX.

# What If the Cluster Has a Fast Flash Tier?

Newer facilities — particularly those built with artificial intelligence workloads in mind — increasingly put an all-flash platform such as VAST, WEKA, or Pure FlashBlade in front of, or instead of, Lustre. This changes things for an architectural reason, not merely because flash is quick: these systems have no separate metadata server. Metadata lives in the same low-latency flash pool as the data, distributed across the namespace, so metadata operations scale with the system rather than funnelling through an MDS pair.

That one design difference is what kills the many-small-files problem. The practical consequence is that much of the advice above inverts:

$TMPDIR still wins in four situations, so do not overcorrect: intermediates you rewrite many times; genuinely latency-bound loops where even a Network File System (NFS) Round-Trip Time (RTT) of tens of microseconds hurts; a shared tier congested by other tenants; and anything requiring Portable Operating System Interface (POSIX) write locking.

What does not change, and this matters:

• Memory per core, walltime, queue behaviour, and scheduler interaction are entirely unaffected.
• MPI halo exchange is a network problem. No storage tier fixes it.
• Hundreds of ranks writing to one GeoPackage still fails. That is file-locking semantics, and over NFS it is frequently worse than on Lustre. Fast storage does not make concurrent writers to a single file safe.
• Per-client throughput is capped by the client’s network card. A node on 100 GbE tops out near 11 GB/s no matter how fast the array is. The levers are nconnect=, NFS over Remote Direct Memory Access (RDMA), or multiple mount points — check what you have with nfsstat -m.
• Flash is expensive, therefore the tier is small, therefore it is usually purged on a 14- or 30-day policy. It is not where results live.

The strongest practical payoff is the object interface. If the platform exposes an Amazon Simple Storage Service (S3) compatible endpoint, then /vsis3/, Zarr, and stackstac-style access work unchanged on-premises:

export AWS_S3_ENDPOINT=vast.internal.example.edu
export AWS_VIRTUAL_HOSTING=FALSE
export AWS_HTTPS=NO
export GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR
gdalinfo /vsis3/elevation/national_1m_dem.tif

One codebase then runs identically on the cluster and in the cloud, which is worth more than any single performance number in this article.

# Recall Is the Only Expensive Transition

Where tape sits behind the disk tiers under Hierarchical Storage Management (HSM), there is exactly one anti-pattern that matters:

The recall storm. A 500-task job array in which each task opens a tape-resident file triggers 500 independent recalls. Cartridges are mounted, sought, and evicted in whatever order requests arrive. A job that should take twenty minutes occupies the archive for a day and stalls every other user on the system.

The fix is to make recall an explicit, separate, batched step, and to have the compute job depend on it:

# 1. stage.sbatch — runs on a data-mover partition, recalls everything at once
#SBATCH --job-name=stage-dem
#SBATCH --partition=copyq
#SBATCH --time=04:00:00
#SBATCH --ntasks=1
sort tile_list.txt > sorted_tiles.txt        # group by cartridge where possible
xargs -a sorted_tiles.txt -P 4 -I{} dmget {} # or the site's stage command

# 2. Submit compute conditional on the stage job succeeding
STAGE_ID=$(sbatch --parsable stage.sbatch)
sbatch --dependency=afterok:$STAGE_ID compute.sbatch

Check residency before assuming — lfs hsm_state on Lustre HSM, dmls/dmattr under Data Migration Facility (DMF), site-specific tooling for High Performance Storage System (HPSS) or Versity. Where the scheduler integrates with staging directly, use it: Slurm’s Burst Buffer (BB) directives and Cray DataWarp (DW) #DW stage_in do this properly, and --tmp= requests node-local scratch capacity.

# How to Tell What You Have

Site documentation is usually vague. Determine it yourself:

stat -f /scratch/ab1234              # filesystem type
mount | grep -E 'lustre|gpfs|nfs|weka|beegfs'
lfs df -h 2>/dev/null                # Lustre present?
nfsstat -m | grep -E 'nconnect|proto' # NFS mount options

Then run the benchmark that actually predicts your workload — not fio, but a geospatial metadata test:

import time, glob, rasterio

paths = sorted(glob.glob("/scratch/ab1234/dem_tiles/*.tif"))[:5000]

t = time.perf_counter()
for p in paths:
    os.stat(p)
print(f"stat  5000 files: {time.perf_counter() - t:.1f}s")

t = time.perf_counter()
for p in paths:
    with rasterio.open(p) as src:
        _ = src.profile          # header only, no pixels
print(f"open  5000 COG headers: {time.perf_counter() - t:.1f}s")

On a metadata-bottlenecked filesystem under load these two numbers diverge badly and the second may run into minutes. On flash both are fast and roughly proportional. Two numbers, five minutes, and you know which half of this section applies to you.

Clusters are frequently hybrid — Lustre for /scratch, a flash tier for an AI partition, tape behind both. The answer is per-path, not per-cluster.

# Interacting With the Scheduler

Now the walkthrough. The workload: zonal statistics over 5,000 one-degree elevation tiles, producing one CSV of per-catchment terrain metrics.

# Step 0 — The Serial Script

Get this working on your laptop against three tiles first. Nothing about the cluster makes debugging easier.

# zonal.py
import sys, rasterio, geopandas as gpd
from rasterstats import zonal_stats

def process_tile(tile_path, catchments_path, out_csv):
    catchments = gpd.read_file(catchments_path)
    with rasterio.open(tile_path) as src:
        subset = catchments[catchments.intersects(box(*src.bounds))]
        if subset.empty:
            return
        stats = zonal_stats(subset, tile_path, stats=["mean", "min", "max", "std"])
    subset.assign(**{k: [s[k] for s in stats] for k in stats[0]}) \
          .drop(columns="geometry") \
          .to_csv(out_csv, index=False)

if __name__ == "__main__":
    process_tile(sys.argv[1], sys.argv[2], sys.argv[3])

# Step 1 — Prove the Environment on a Compute Node

Do not debug inside a batch job. Get an interactive shell on a real compute node and run exactly one tile:

srun --pty --time=00:30:00 --mem=8G --cpus-per-task=1 bash
module load gdal/3.8.0
apptainer exec --bind /scratch:/scratch geo.sif \
  python zonal.py /scratch/ab1234/dem_tiles/tile_0001.tif catchments.gpkg /tmp/t1.csv

If this works, everything that follows is mechanical. If it does not, nothing that follows will help.

# Step 2 — Your First Batch Job

#!/bin/bash
#SBATCH --job-name=zonal-one        # appears in squeue
#SBATCH --account=xy00              # which allocation to charge
#SBATCH --partition=normal
#SBATCH --time=00:20:00             # walltime; job is KILLED at this point
#SBATCH --ntasks=1                  # one process
#SBATCH --cpus-per-task=1           # one core for it
#SBATCH --mem=8G                    # memory for the whole job
#SBATCH --output=logs/%x_%j.out     # %x = job name, %j = job id
#SBATCH --error=logs/%x_%j.err

module load gdal/3.8.0
export OMP_NUM_THREADS=1            # do not let libraries fight each other

srun apptainer exec --bind /scratch:/scratch geo.sif \
  python zonal.py "$TILE" catchments.gpkg "out/$(basename $TILE .tif).csv"

sbatch --export=TILE=/scratch/ab1234/dem_tiles/tile_0001.tif zonal_one.sbatch
squeue -u $USER

When it finishes, ask what it actually used:

sacct -j 4821993 --format=JobID,JobName,Elapsed,TotalCPU,MaxRSS,State
#   JobID   JobName    Elapsed   TotalCPU     MaxRSS      State
# 4821993   zonal-one  00:00:47  00:00:44      1.9G   COMPLETED

MaxRSS of 1.9 GB against a request of 8 GB means the request was three times too large. This matters more than it looks: the scheduler cannot start your job until a node has your requested memory free, so over-requesting is a direct cause of long queue waits. Right-size from measurement, add perhaps 30% headroom, and your throughput improves without any code change.

# Step 3 — Job Arrays, Which Solve Most Spatial Problems

A job array submits one script N times with a different index. This is the answer for the large majority of geospatial batch work.

#!/bin/bash
#SBATCH --job-name=zonal-array
#SBATCH --account=xy00
#SBATCH --time=00:30:00
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=1
#SBATCH --mem=3G
#SBATCH --array=0-499%50            # 500 tasks, at most 50 running at once
#SBATCH --output=logs/%x_%A_%a.out   # %A = array job id, %a = task index

module load gdal/3.8.0
export OMP_NUM_THREADS=1

srun apptainer exec --bind /scratch:/scratch geo.sif \
  python zonal_chunk.py --index $SLURM_ARRAY_TASK_ID --nchunks 500

The script slices its own share of the work from the index — no coordination, no communication:

# zonal_chunk.py
import argparse, glob, os
from zonal import process_tile

ap = argparse.ArgumentParser()
ap.add_argument("--index", type=int, required=True)
ap.add_argument("--nchunks", type=int, required=True)
args = ap.parse_args()

tiles = sorted(glob.glob("/scratch/ab1234/dem_tiles/*.tif"))
mine = tiles[args.index::args.nchunks]     # strided: even work per task
print(f"task {args.index}: {len(mine)} tiles", flush=True)

for t in mine:
    out = f"out/{os.path.basename(t)[:-4]}.csv"
    if os.path.exists(out):                # idempotent: safe to re-run
        continue
    process_tile(t, "catchments.gpkg", out)

Two details earn their keep. The %50 throttle limits how many tasks run concurrently, which keeps you from monopolising a partition and, on a Lustre system, from creating your own metadata storm. And the os.path.exists check makes the job idempotent, so recovering from a partial failure is sbatch --array=17,203,488 zonal_array.sbatch rather than a full re-run.

# Step 4 — Using a Whole Node

Where each task is chunky enough to justify it, ask for multiple cores and use them:

#SBATCH --cpus-per-task=16
#SBATCH --mem=48G

from joblib import Parallel, delayed
import os

n = int(os.environ.get("SLURM_CPUS_PER_TASK", 1))
Parallel(n_jobs=n, backend="loky")(
    delayed(process_tile)(t, "catchments.gpkg", f"out/{os.path.basename(t)[:-4]}.csv")
    for t in mine
)

The trap here is thread oversubscription. GDAL, NumPy’s Basic Linear Algebra Subprograms (BLAS) backend, and Open Multi-Processing (OpenMP) will each happily spawn 16 threads inside each of your 16 workers, producing 256 threads fighting over 16 cores. Pin them:

export OMP_NUM_THREADS=1
export OPENBLAS_NUM_THREADS=1
export MKL_NUM_THREADS=1
export GDAL_NUM_THREADS=1
export GDAL_CACHEMAX=512            # MB, per process — not per node

# Step 5 — Checkpointing Against Walltime

Long jobs get killed. Catch the warning signal and save state:

#SBATCH --time=04:00:00
#SBATCH --signal=B:USR1@300         # send USR1 to the batch step 5 min before the kill

import signal, json

state = {"done": []}

def checkpoint(signum, frame):
    with open("checkpoint.json", "w") as f:
        json.dump(state, f)
    raise SystemExit(0)

signal.signal(signal.SIGUSR1, checkpoint)

Combined with the idempotence check from Step 3, this makes a requeued job resume rather than restart.

# Distributing Across Nodes

Two models, and the choice is determined by whether your tasks need to talk to each other.

# Loosely Coupled: dask-jobqueue

dask-jobqueue submits Slurm jobs as Dask workers, so you keep writing Dask code and the cluster grows and shrinks underneath it:

from dask_jobqueue import SLURMCluster
from dask.distributed import Client
import dask.dataframe as dd

cluster = SLURMCluster(
    queue="normal",
    account="xy00",
    cores=16,                  # cores per Slurm job
    processes=8,               # split into 8 workers of 2 cores each
    memory="48GB",
    walltime="02:00:00",
    local_directory="$TMPDIR", # spill to node-local NVMe, not shared storage
    job_extra_directives=["--export=ALL", "--output=logs/dask-%j.out"],
    job_script_prologue=[
        "module load gdal/3.8.0",
        "export OMP_NUM_THREADS=1",
    ],
)
cluster.adapt(minimum_jobs=1, maximum_jobs=25)   # up to 400 cores on demand
client = Client(cluster)

ddf = dd.read_parquet("/scratch/ab1234/catchments_partitioned/")
result = ddf.groupby("catchment_id").elevation.mean().compute()

processes=8 with cores=16 is deliberate: eight separate Python processes sidestep the Global Interpreter Lock (GIL) that would otherwise limit a single 16-thread worker, and local_directory pointed at $TMPDIR keeps Dask’s spill traffic off the shared filesystem.

On clusters where compute nodes cannot open arbitrary network connections, adapt() will fail because workers cannot reach the scheduler. The robust alternative is to run the whole Dask cluster inside one allocation, coordinating through a file on shared storage:

#!/bin/bash
#SBATCH --nodes=8 --exclusive --time=02:00:00 --account=xy00

SCHED=/scratch/ab1234/sched-$SLURM_JOB_ID.json
srun -N1 -n1 dask scheduler --scheduler-file $SCHED &
sleep 15
srun -N7 -n56 dask worker --scheduler-file $SCHED \
  --nthreads 2 --memory-limit 5GB --local-directory $TMPDIR &
sleep 20
python analysis.py --scheduler-file $SCHED

# Tightly Coupled: mpi4py and Halo Exchange

Some spatial operations genuinely need neighbours. Slope, aspect, flow direction, flow accumulation, and any focal window operation compute each output cell from the cells around it. Split a Digital Elevation Model (DEM) into row blocks across ranks and every block needs the last row of the block above and the first row of the block below — the halo, or ghost rows.

This is what Dask handles clumsily and MPI handles well:

from mpi4py import MPI
import numpy as np, rasterio
from scipy.ndimage import sobel

comm = MPI.COMM_WORLD
rank, size = comm.Get_rank(), comm.Get_size()
HALO = 1

with rasterio.open("/scratch/ab1234/national_1m_dem.tif") as src:
    nrows, ncols = src.height, src.width
    rows_each = nrows // size
    r0 = rank * rows_each
    r1 = nrows if rank == size - 1 else r0 + rows_each

    # Read own block plus a halo row on each interior side
    read_start = max(0, r0 - HALO)
    read_stop = min(nrows, r1 + HALO)
    block = src.read(1, window=((read_start, read_stop), (0, ncols))).astype("f4")
    res = src.res[0]

up_ghost = np.empty(ncols, dtype="f4")
down_ghost = np.empty(ncols, dtype="f4")

# Exchange boundaries: send my top row up, receive my neighbour's bottom row
comm.Sendrecv(sendbuf=np.ascontiguousarray(block[HALO]),      dest=rank - 1 if rank > 0 else MPI.PROC_NULL,
              recvbuf=down_ghost,                             source=rank + 1 if rank < size - 1 else MPI.PROC_NULL)
comm.Sendrecv(sendbuf=np.ascontiguousarray(block[-HALO - 1]), dest=rank + 1 if rank < size - 1 else MPI.PROC_NULL,
              recvbuf=up_ghost,                               source=rank - 1 if rank > 0 else MPI.PROC_NULL)

# Compute slope on the padded block, then discard the halo rows
dz_dx = sobel(block, axis=1) / (8 * res)
dz_dy = sobel(block, axis=0) / (8 * res)
slope = np.degrees(np.arctan(np.hypot(dz_dx, dz_dy)))
interior = slope[HALO:-HALO] if 0 < rank < size - 1 else slope

# Each rank writes its own file — never a shared handle
np.save(f"/scratch/ab1234/slope_rank{rank:04d}.npy", interior)
comm.Barrier()
if rank == 0:
    print("all ranks complete", flush=True)

Launched with:

#SBATCH --nodes=8
#SBATCH --ntasks-per-node=48
#SBATCH --exclusive
module load openmpi/4.1.5 gdal/3.8.0
srun --mpi=pmix python slope_mpi.py

Note the last two rules again, because they are the ones people break: one output file per rank, and a barrier before any rank assumes the others are finished. Writing into a single shared GeoTIFF from 384 ranks will corrupt it or hang, on any filesystem, at any speed.

# Proving It Was Worth It

Measure two things. Strong scaling holds the problem fixed and adds cores — ideal is linear, reality tails off as serial fractions and communication dominate. Weak scaling grows the problem with the cores, which is the more honest measure for tiled geospatial work.

seff 4821993
# Job Wall-clock time: 00:41:12
# CPU Efficiency: 27.3% of 05:29:36 core-walltime
# Memory Utilized: 2.11 GB (Est. 8.00 GB)

Twenty-seven percent CPU efficiency means roughly three-quarters of the core-hours you were charged did nothing. The arithmetic is unforgiving: 400 ranks at 25% efficiency burns four times the allocation of 100 ranks at 90% and finishes at a similar wall-clock time. Amdahl’s law sets the ceiling — if 5% of your pipeline is inherently serial, no core count takes you past a 20× speedup.

# Pitfalls, Collected

# Your First Week, In Order

1. Set up ~/.ssh/config with ControlMaster. Run sinfo -s and module avail gdal.
2. Run the two-number storage benchmark from earlier. Decide which half of the storage section applies to you.
3. Build an Apptainer image with your geospatial stack. Verify it on an interactive node with one tile.
4. Submit one batch job. Read sacct output. Right-size memory and walltime.
5. Convert it to a job array with a %N throttle and an idempotence check.
6. Only if the work is genuinely coupled, reach for MPI. Most geospatial work never needs it.
7. Run seff on a completed array. If efficiency is under 50%, fix that before asking for more cores.

The step most people skip is the second, and it is the one that determines whether the rest of the advice is right or exactly backwards.

Related reading: Dask-GeoPandas: Parallel Spatial Processing Across CPU Cores · Rasterio Windowed Reading: Processing 100 GB Rasters on 8 GB of RAM · cuSpatial: GPU-Accelerated Spatial Analytics with RAPIDS