Short answer

A geological map database in PostGIS is not a single layer of coloured polygons. It is a small relational model: a polygon table of map units, a linework table of contacts and faults, a point table of structural measurements, and a lithostratigraphic lookup table that every spatial feature references by a stable code. PostGIS adds the geometry types, spatial indexing, on-server validation, and SQL that let several analysts, several map sheets, and a web service all read the same authoritative data without duplicating it. The work is mostly schema design and topology discipline; the spatial functions are the easy part.

Why a database, not a folder of shapefiles

Most geological mapping starts as digitised polygons in QGIS, exported as shapefiles or a GeoPackage per sheet. That is fine for one mapper on one sheet. It breaks the moment you have ten 1:50,000 sheets, three field geologists editing concurrently, and a requirement to symbolise the whole region by the same legend.

PostGIS solves four specific problems a file collection cannot:

  • Concurrent editing. PostgreSQL handles row-level locking, so two geologists can edit different sheets, or different features of the same sheet, at once. A shapefile or GeoPackage is effectively single-writer.
  • Referential integrity. A FOREIGN KEY from your polygon table to a strat_units lexicon means you physically cannot save a polygon tagged with a unit code that does not exist. No more "Jurassic" in one sheet and "jurassic" in the next.
  • Server-side logic. Validity, area, adjacency, and dissolve operations run as SQL against the live data, not as a manual QGIS step someone forgets.
  • One source, many products. The same tables feed a print map, a vector-tile web service, and an attribute query, with no export-and-drift.

A schema that works

A defensible minimum for a geological map is four tables. Use geometry(MultiPolygon, 25832) rather than generic geometry so the column enforces type and SRID; here 25832 is ETRS89 / UTM zone 32N (EPSG:25832), a typical projected CRS for central Europe — pick the one that matches your area of work.

CREATE TABLE strat_units (
    unit_code   text PRIMARY KEY,        -- e.g. 'J2_lst'
    unit_name   text NOT NULL,           -- 'Middle Jurassic limestone'
    lithology   text,                    -- 'limestone'
    age_period  text,                    -- 'Jurassic'
    age_min_ma  numeric,                 -- 163.5
    age_max_ma  numeric,                 -- 174.1
    map_color   text                     -- hex, optional cartographic hint
);

CREATE TABLE map_units (
    gid         serial PRIMARY KEY,
    unit_code   text NOT NULL REFERENCES strat_units(unit_code),
    sheet       text,                     -- source map sheet id
    geom        geometry(MultiPolygon, 25832) NOT NULL
);

CREATE TABLE contacts (
    gid         serial PRIMARY KEY,
    contact_type text,                    -- 'depositional','fault','unconformity'
    confidence  text,                     -- 'observed','inferred','concealed'
    geom        geometry(MultiLineString, 25832) NOT NULL
);

CREATE TABLE structure_pts (
    gid         serial PRIMARY KEY,
    obs_type    text,                     -- 'bedding','foliation','joint'
    dip_dir     smallint,                 -- 0-360
    dip         smallint,                 -- 0-90
    geom        geometry(Point, 25832) NOT NULL
);

Two design choices carry most of the value. First, contacts are their own table, not polygon boundaries. A geological contact between two units is a real object with its own attributes (observed vs inferred, fault vs depositional), and you want to edit and symbolise it independently of the polygons it happens to separate. Second, everything joins back to strat_units by code, so age, lithology and colour are defined once. To draw the map you join: SELECT m.geom, s.map_color, s.age_period FROM map_units m JOIN strat_units s USING (unit_code);.

Loading data

Import existing sheets with GDAL's ogr2ogr, which writes straight into PostGIS and reprojects on the way in:

ogr2ogr -f PostgreSQL \
  PG:"dbname=geology user=mapper" \
  sheet_7421.gpkg map_units \
  -nln map_units -append \
  -t_srs EPSG:25832 \
  -nlt PROMOTE_TO_MULTI \
  -lco GEOMETRY_NAME=geom -lco FID=gid

-nlt PROMOTE_TO_MULTI avoids the common failure where a layer holds both Polygon and MultiPolygon and the typed column rejects the singles. -t_srs reprojects; do not confuse this with merely declaring an SRID. If your source is genuinely already in EPSG:25832 and only the metadata is missing, use -a_srs to assign instead. After every bulk load, run ANALYZE map_units; so the query planner has fresh statistics.

Indexing and performance

A geometry column with no spatial index forces a sequential scan on every ST_Intersects or bounding-box filter. Create a GiST index on each:

CREATE INDEX map_units_geom_gix ON map_units USING GIST (geom);
CREATE INDEX contacts_geom_gix  ON contacts  USING GIST (geom);
CREATE INDEX map_units_unit_idx ON map_units (unit_code);   -- B-tree for the join

The GiST index accelerates the bounding-box stage of spatial predicates; PostGIS then refines with the exact geometry test. For a regional map of a few hundred thousand polygons this is the difference between sub-second map rendering and timeouts. Confirm the planner is actually using the index with EXPLAIN ANALYZE on a representative tile query.

Topology: the part that bites

Geological maps must be a clean planar partition — no gaps, no overlaps between adjacent units. PostGIS gives you two ways to enforce this.

The lightweight approach is validation queries you run after editing:

-- Invalid geometries (self-intersections from sloppy digitising)
SELECT gid, ST_IsValidReason(geom) FROM map_units WHERE NOT ST_IsValid(geom);

-- Overlaps between units (should return nothing)
SELECT a.gid, b.gid
FROM map_units a JOIN map_units b
  ON a.gid < b.gid AND ST_Overlaps(a.geom, b.geom);

Repair self-intersections with UPDATE map_units SET geom = ST_MakeValid(geom) WHERE NOT ST_IsValid(geom);. The heavier, more correct approach is the postgis_topology extension, which stores edges and faces explicitly so a shared boundary exists once and gaps are structurally impossible. It costs more setup but pays off when many sheets must tile seamlessly.

Sheet edges are where most problems hide. When two 1:50,000 sheets were digitised separately, their shared boundary rarely matches to the millimetre, leaving sliver polygons and hairline gaps. Snap them before merging:

UPDATE map_units a
SET geom = ST_Snap(a.geom, b.geom, 0.5)   -- 0.5 m tolerance, projected CRS
FROM map_units b
WHERE a.sheet <> b.sheet AND ST_DWithin(a.geom, b.geom, 0.5);

Choose the tolerance from your mapping accuracy, not arbitrarily: at 1:50,000, half a metre on the ground is far below line weight, so 0.5–1 m is safe; a tolerance larger than your thinnest real polygon will eat genuine features.

Validation before you trust the database

  • Coverage: SELECT ST_Area(ST_Union(geom)) FROM map_units; should equal the mapped extent. A shortfall means gaps.
  • Orphan codes: the foreign key prevents these, but verify after a raw COPY load that bypassed it.
  • Age sanity: age_min_ma < age_max_ma for every unit; flag rows that violate it.
  • Contact–polygon agreement: every fault in contacts should lie on or near a polygon boundary; ST_DWithin against unioned boundaries surfaces strays.

Bathyl perspective

We treat a geological map database as a model of stratigraphic reality, not a styled picture. The schema — typed geometry columns, a single lexicon, contacts as first-class objects, enforced validity — is what lets a map survive new sheets, a second geologist, and a later web product without quietly drifting out of sync. Getting the tables right is cheaper than fixing a region-wide overlap a year later.

Related reading

Sources