Short answer

Command line GIS — chiefly GDAL/OGR and PostGIS — earns its place the moment an operation has to be repeated, scaled, or audited. A single reprojection is faster in QGIS; reprojecting 400 lidar tiles, re-running the same cleanup every time a client sends a fresh export, or proving exactly what was done to a deliverable is where the terminal wins decisively. The payoff is not speed for its own sake: it is that a command is rerunnable, inspectable, and version-controllable, so the work survives new data and a new analyst. Below are the recipes a consultant actually reaches for.

The decision: when to leave the GUI

Stay in QGIS/ArcGIS when the task is exploratory, visual, or one-off. Move to the command line when any of these is true:

  • Volume: the same operation across many files (tiles, monthly exports, multiple clients).
  • Repetition: the operation will recur as source data updates.
  • Auditability: you must be able to show, exactly, what transformed the data.
  • Scale: the dataset is too large to load comfortably, or the work belongs in a database.
  • Integration: the step is part of an automated pipeline or scheduled job.

The cost is upfront care: command line tools do exactly what you say, including the wrong thing, silently. So the workflow is always inspect → prototype on one file → run the batch → validate.

Inspect before you touch anything

Never convert blind. Read the source first.

Raster:

gdalinfo dem.tif

Tells you the CRS, corner coordinates, pixel size, band count, data type, NoData value, and whether overviews exist. The NoData value matters: warp or resample without it set correctly and the fill value bleeds into your data.

Vector:

ogrinfo -al -so data.gpkg

-so (summary only) prints the layer SRS, geometry type, extent, feature count, and the full field list with types — without dumping every feature. This is how you catch a missing CRS or a schema change before it propagates.

Just the CRS, machine-readable:

gdalsrsinfo -o epsg data.gpkg

Reprojection, done right

Raster reprojection with gdalwarp. Resampling choice is not cosmetic: use a continuous method for elevation/imagery and nearest-neighbour for categorical data, or you will invent values that never existed.

# DEM (continuous): bilinear or cubic
gdalwarp -s_srs EPSG:4326 -t_srs EPSG:32631 -r bilinear \
  -co COMPRESS=DEFLATE -co TILED=YES dem_4326.tif dem_utm.tif

# Land-cover classes (categorical): nearest, never bilinear
gdalwarp -t_srs EPSG:32631 -r near landcover.tif landcover_utm.tif

Vector reprojection with ogr2ogr — note the distinction that bites people: -t_srs reprojects (recomputes coordinates), -a_srs only assigns (relabels, no movement):

# Reproject from a known source CRS
ogr2ogr -s_srs EPSG:4326 -t_srs EPSG:32631 out_utm.gpkg in.gpkg

# Assign a missing CRS without moving the data
ogr2ogr -a_srs EPSG:32631 fixed.gpkg unknown_crs.shp

When datums differ (NAD27↔NAD83, OSGB36↔WGS84), check that PROJ picked the right transformation, because a wrong datum operation produces confidently wrong coordinates: projinfo -s EPSG:27700 -t EPSG:4326 --summary.

A worked batch: 400 DEM tiles to one analytical mosaic

A common consulting task: a client ships hundreds of DEM tiles in WGS84, and you need one reprojected, compressed, cloud-friendly DEM in UTM for slope analysis.

# 1. Reproject each tile (continuous data -> bilinear)
mkdir reproj
for f in tiles/*.tif; do
  gdalwarp -t_srs EPSG:32631 -r bilinear -co COMPRESS=DEFLATE \
    "$f" "reproj/$(basename "$f")"
done

# 2. Build a virtual mosaic (no data duplication)
gdalbuildvrt mosaic.vrt reproj/*.tif

# 3. Materialize a Cloud Optimized GeoTIFF with overviews
gdal_translate mosaic.vrt dem_mosaic_cog.tif -of COG \
  -co COMPRESS=DEFLATE -co OVERVIEW_RESAMPLING=BILINEAR

# 4. Derive slope directly from the mosaic
gdaldem slope dem_mosaic_cog.tif slope.tif -compute_edges

gdalbuildvrt is the key efficiency trick: the VRT is a tiny XML index over the tiles, so you mosaic terabytes without copying them, then materialize only the final product.

PostGIS from the terminal

When several people or apps need to query the same spatial truth, files stop scaling and a database earns its keep. Load with ogr2ogr:

ogr2ogr -f PostgreSQL "PG:dbname=gis user=me" survey.gpkg \
  -t_srs EPSG:32631 \
  -lco GEOMETRY_NAME=geom -lco FID=id -nlt PROMOTE_TO_MULTI -nln survey

Then index and query in SQL. A spatial index is mandatory before any join over real data:

CREATE INDEX survey_geom_idx ON survey USING GIST (geom);

-- Find samples within 500 m of a fault, in projected meters
SELECT s.id
FROM survey s
JOIN faults f ON ST_DWithin(s.geom, f.geom, 500);

Two PostGIS functions consultants confuse constantly: ST_SetSRID(geom, 32631) only stamps the SRID (assign), while ST_Transform(geom, 4326) reprojects (recompute). ST_DWithin is preferred over ST_Distance(...) < x because it uses the spatial index.

Common pitfalls and why they happen

  • Resampling categorical rasters with bilinear. Averaging class codes (e.g., averaging "granite=3" and "shale=5" into "4") invents nonexistent classes. Use -r near for thematic data.
  • ST_SetSRID when transformation was needed. It relabels without moving, so the geometry stays in the old position under a new label. Use ST_Transform when the coordinates must change.
  • Assuming success from an output file existing. A truncated or empty output still creates a file. Validate counts and extents.
  • Forgetting NoData on warp. The fill value leaks into the data edges. Confirm NoData in gdalinfo and pass -dstnodata when needed.
  • Batch converting through changing schemas. If tile 200 has a different field set, a naive loop may drop columns silently. Inspect a sample of inputs, not just the first.
  • No spatial index in PostGIS. A spatial join without a GIST index does a full cross product and crawls.

QA and validation

  • Re-run gdalinfo/ogrinfo on outputs and compare CRS, extent, band/feature counts, and NoData against the source.
  • For vectors, compare feature counts: ogrinfo -so before and after; a drop means dropped or invalid geometries.
  • Round-trip a known coordinate through any reprojection; it should return to within sub-millimeter.
  • Check geometry validity (SELECT count(*) FROM t WHERE NOT ST_IsValid(geom); should be 0).
  • Keep the commands themselves in a script under version control — they are the documentation of what was done.

Bathyl perspective

We move work to the command line when GIS becomes infrastructure rather than a one-off map. The goal is never cleverness for its own sake; it is a pipeline that can be rerun on next quarter's data by a different analyst and produce the same, defensible result. A script under version control is both the work and its audit trail — which is exactly what a consulting deliverable should be able to stand on.

Related reading

Sources