Geospatial Data Processing with PostGIS and GeoServer

The Need for Spatial Data Infrastructure

Most data pipelines treat location as a flat attribute -- a latitude/longitude pair stored in two float columns. This works until you need to answer questions that are inherently spatial: which facilities are within a 30-minute walk? Where do service coverage areas overlap? Which parcels fall inside a flood zone? At that point, you need a proper spatial data infrastructure that can index, query, and serve geographic data at scale.

PostGIS and GeoServer form the backbone of open-source spatial infrastructure used by governments, utilities, and logistics companies worldwide. PostGIS extends PostgreSQL with geometry types, spatial indexes, and hundreds of analytical functions. GeoServer sits on top, exposing those datasets as standardised OGC web services -- WMS for rendered map tiles, WFS for raw feature data, and WMTS for cached tile pyramids.

This article walks through the full pipeline: setting up PostGIS with proper indexing, ingesting spatial data at scale, writing advanced spatial queries, configuring GeoServer for tile serving, and integrating with MapLibre GL on the front end. The final sections cover a real-world healthcare isochrone analysis for Hamburg and production operations patterns.

PostGIS Foundation

PostGIS is not just a set of functions bolted onto PostgreSQL -- it introduces new data types, new index types, and a catalogue of spatial reference systems that the query planner understands natively. Getting the foundation right matters because spatial queries that hit millions of geometries are extremely sensitive to type choices and index strategies.

Extension Setup

-- Enable PostGIS and related extensions
CREATE EXTENSION IF NOT EXISTS postgis;
CREATE EXTENSION IF NOT EXISTS postgis_topology;
CREATE EXTENSION IF NOT EXISTS postgis_raster;
CREATE EXTENSION IF NOT EXISTS fuzzystrmatch;
CREATE EXTENSION IF NOT EXISTS postgis_tiger_geocoder;

-- Verify installation
SELECT PostGIS_Full_Version();
-- POSTGIS="3.4.2" PGSQL="160" GEOS="3.12.1" PROJ="9.3.1" LIBXML="2.11.6"

Spatial Reference Systems (SRID)

Every geometry in PostGIS is associated with a Spatial Reference System Identifier (SRID). The two most common are EPSG:4326 (WGS 84 -- longitude/latitude in degrees) and EPSG:3857 (Web Mercator -- meters, used by most web maps). For regional analysis in Germany, EPSG:25832 (ETRS89 / UTM zone 32N) provides metric precision without global distortion.

-- Check available SRIDs
SELECT srid, auth_name, auth_srid, proj4text
FROM spatial_ref_sys
WHERE auth_name = 'EPSG' AND auth_srid IN (4326, 3857, 25832);

-- Transform between coordinate systems
SELECT ST_Transform(
  ST_SetSRID(ST_MakePoint(9.9937, 53.5511), 4326),
  25832
) AS hamburg_utm;

Geometry vs Geography Types

-- Create a table with geometry type (projected coordinates)
CREATE TABLE facilities (
  id          SERIAL PRIMARY KEY,
  name        TEXT NOT NULL,
  category    TEXT NOT NULL,
  geom        geometry(Point, 25832) NOT NULL,
  created_at  TIMESTAMPTZ DEFAULT now()
);

-- Create a table with geography type (lat/lon)
CREATE TABLE global_sensors (
  id          SERIAL PRIMARY KEY,
  sensor_id   TEXT UNIQUE NOT NULL,
  location    geography(Point, 4326) NOT NULL,
  last_seen   TIMESTAMPTZ DEFAULT now()
);

Indexing Strategies

Spatial indexes are the single most important factor for query performance. PostGIS supports GiST (Generalized Search Tree) and SP-GiST indexes. GiST is the default and handles bounding-box overlap queries efficiently. SP-GiST uses space-partitioning trees and can be faster for point-only datasets with highly skewed distributions.

-- GiST index (default, general purpose)
CREATE INDEX idx_facilities_geom
  ON facilities USING GIST (geom);

-- SP-GiST index (better for point clouds)
CREATE INDEX idx_sensors_location
  ON global_sensors USING SPGIST (location);

-- Partial index for active facilities only
CREATE INDEX idx_active_facilities_geom
  ON facilities USING GIST (geom)
  WHERE category != 'decommissioned';

-- Verify index usage
EXPLAIN (ANALYZE, BUFFERS)
SELECT id, name
FROM facilities
WHERE ST_DWithin(
  geom,
  ST_Transform(ST_SetSRID(ST_MakePoint(9.9937, 53.5511), 4326), 25832),
  5000
);

Data Ingestion Pipelines

Spatial data arrives in dozens of formats -- Shapefiles, GeoJSON, GeoPackage, KML, CSV with coordinates, and raw GPS traces. A robust ingestion pipeline must handle format conversion, coordinate transformation, validation, and bulk loading without choking on multi-gigabyte datasets.

Loading Shapefiles with ogr2ogr

# Load a shapefile into PostGIS with coordinate transformation
ogr2ogr -f "PostgreSQL" \
  PG:"host=localhost dbname=geodata user=geouser password=secret" \
  /data/shapefiles/hamburg_buildings.shp \
  -nln public.buildings \
  -nlt PROMOTE_TO_MULTI \
  -t_srs EPSG:25832 \
  -lco GEOMETRY_NAME=geom \
  -lco FID=id \
  -lco PRECISION=NO \
  --config PG_USE_COPY YES \
  -progress

# Load GeoJSON with append mode
ogr2ogr -f "PostgreSQL" \
  PG:"host=localhost dbname=geodata user=geouser" \
  /data/geojson/points_of_interest.geojson \
  -nln public.poi \
  -append \
  -t_srs EPSG:25832

# Load GeoPackage selecting specific layers
ogr2ogr -f "PostgreSQL" \
  PG:"host=localhost dbname=geodata user=geouser" \
  /data/packages/transport.gpkg \
  -sql "SELECT * FROM bus_stops WHERE status = 'active'" \
  -nln public.bus_stops

Python Ingestion with GeoPandas and Fiona

import geopandas as gpd
import fiona
from sqlalchemy import create_engine

engine = create_engine(
    "postgresql://geouser:secret@localhost:5432/geodata"
)

# Read shapefile with GeoPandas
gdf = gpd.read_file("/data/shapefiles/hamburg_districts.shp")

# Inspect and reproject
print(f"Original CRS: {gdf.crs}")
print(f"Records: {len(gdf)}")
gdf = gdf.to_crs(epsg=25832)

# Validate geometries
invalid = gdf[~gdf.geometry.is_valid]
print(f"Invalid geometries: {len(invalid)}")
gdf.geometry = gdf.geometry.buffer(0)  # fix common issues

# Write to PostGIS
gdf.to_postgis(
    name="districts",
    con=engine,
    schema="public",
    if_exists="replace",
    index=True,
    index_label="id",
    dtype={"geometry": "geometry"},
)

print("Districts loaded successfully.")

Bulk Loading with Fiona for Large Datasets

import fiona
from shapely.geometry import shape, mapping
from shapely.ops import transform
import pyproj
import psycopg2

# Define coordinate transformation
project = pyproj.Transformer.from_crs(
    "EPSG:4326", "EPSG:25832", always_xy=True
).transform

conn = psycopg2.connect(
    host="localhost", dbname="geodata",
    user="geouser", password="secret"
)
cur = conn.cursor()

batch = []
batch_size = 5000

with fiona.open("/data/large_dataset.gpkg", layer="parcels") as src:
    for i, feature in enumerate(src):
        geom = shape(feature["geometry"])
        geom_projected = transform(project, geom)
        props = feature["properties"]

        batch.append((
            props.get("parcel_id"),
            props.get("area_sqm"),
            geom_projected.wkb_hex,
        ))

        if len(batch) >= batch_size:
            cur.executemany(
                """INSERT INTO parcels (parcel_id, area_sqm, geom)
                   VALUES (%s, %s, ST_GeomFromWKB(%s::geometry, 25832))
                   ON CONFLICT (parcel_id) DO UPDATE
                   SET area_sqm = EXCLUDED.area_sqm,
                       geom = EXCLUDED.geom""",
                batch,
            )
            conn.commit()
            batch.clear()
            print(f"Loaded {i + 1} features...")

    if batch:
        cur.executemany(
            """INSERT INTO parcels (parcel_id, area_sqm, geom)
               VALUES (%s, %s, ST_GeomFromWKB(%s::geometry, 25832))
               ON CONFLICT (parcel_id) DO UPDATE
               SET area_sqm = EXCLUDED.area_sqm,
                   geom = EXCLUDED.geom""",
            batch,
        )
        conn.commit()

cur.close()
conn.close()
print("Bulk load complete.")

Spatial Query Patterns

PostGIS provides hundreds of spatial functions. The patterns below cover the most common analytical operations you will need in production -- from simple proximity searches to complex coverage gap detection.

Buffer Analysis

-- Find all residential buildings within 500m of a hospital
SELECT b.id, b.address, b.building_type,
       ST_Distance(b.geom, h.geom) AS distance_m
FROM buildings b
JOIN hospitals h ON ST_DWithin(b.geom, h.geom, 500)
WHERE b.building_type = 'residential'
  AND h.name = 'UKE Hamburg'
ORDER BY distance_m;

-- Create buffer zones around industrial sites
SELECT s.id, s.name,
       ST_Buffer(s.geom, 200) AS buffer_200m,
       ST_Buffer(s.geom, 500) AS buffer_500m
FROM industrial_sites s;

Point-in-Polygon and Spatial Joins

-- Count facilities per district (spatial join)
SELECT d.name AS district,
       d.population,
       COUNT(f.id) AS facility_count,
       ROUND(d.population::numeric / NULLIF(COUNT(f.id), 0)) AS pop_per_facility
FROM districts d
LEFT JOIN facilities f ON ST_Contains(d.geom, f.geom)
GROUP BY d.name, d.population
ORDER BY pop_per_facility DESC NULLS FIRST;

-- Point-in-polygon: which zone does a coordinate fall in?
SELECT z.zone_name, z.zone_type
FROM zones z
WHERE ST_Contains(
  z.geom,
  ST_Transform(ST_SetSRID(ST_MakePoint(9.9937, 53.5511), 4326), 25832)
);

Nearest Neighbor Search

-- Find the 5 nearest pharmacies to a given point
-- Uses the <-> operator for index-assisted KNN
SELECT p.id, p.name, p.address,
       ST_Distance(
         p.geom,
         ST_Transform(ST_SetSRID(ST_MakePoint(10.0, 53.55), 4326), 25832)
       ) AS distance_m
FROM pharmacies p
ORDER BY p.geom <->
  ST_Transform(ST_SetSRID(ST_MakePoint(10.0, 53.55), 4326), 25832)
LIMIT 5;

Coverage Gap Detection

-- Identify areas NOT covered by any facility within 2km
WITH facility_coverage AS (
  SELECT ST_Union(ST_Buffer(f.geom, 2000)) AS covered_area
  FROM facilities f
  WHERE f.category = 'healthcare'
),
city_boundary AS (
  SELECT ST_Union(geom) AS boundary
  FROM districts
)
SELECT ST_Difference(cb.boundary, fc.covered_area) AS gap_area,
       ST_Area(ST_Difference(cb.boundary, fc.covered_area)) / 1e6 AS gap_km2
FROM city_boundary cb, facility_coverage fc;

Isochrone Generation

-- Generate walking isochrones using pgRouting
-- Requires a routable network graph
SELECT id, the_geom, agg_cost
FROM pgr_drivingDistance(
  'SELECT id, source, target, cost_walk AS cost,
          reverse_cost_walk AS reverse_cost
   FROM road_network',
  (SELECT id FROM road_network_vertices_pgr
   ORDER BY the_geom <-> ST_SetSRID(ST_MakePoint(9.9937, 53.5511), 4326)
   LIMIT 1),
  ARRAY[900, 1800, 2700],  -- 15, 30, 45 minutes in seconds
  directed := false
);

-- Build isochrone polygons from reachable nodes
WITH reachable AS (
  SELECT dd.agg_cost, v.the_geom
  FROM pgr_drivingDistance(
    'SELECT id, source, target, cost_walk AS cost,
            reverse_cost_walk AS reverse_cost
     FROM road_network',
    42,  -- source vertex id
    ARRAY[900, 1800, 2700],
    directed := false
  ) dd
  JOIN road_network_vertices_pgr v ON dd.node = v.id
)
SELECT
  CASE
    WHEN agg_cost <= 900 THEN '15 min'
    WHEN agg_cost <= 1800 THEN '30 min'
    ELSE '45 min'
  END AS isochrone_band,
  ST_ConcaveHull(ST_Collect(the_geom), 0.7) AS geom
FROM reachable
GROUP BY 1;

GeoServer Configuration

GeoServer publishes PostGIS layers as OGC-compliant web services. The setup follows a hierarchy: Workspace (logical grouping) contains Stores (database connections) which contain Layers (individual tables or views). Styles are applied via SLD (Styled Layer Descriptor) documents.

Docker Deployment

# docker-compose.yml
version: "3.8"
services:
  postgis:
    image: postgis/postgis:16-3.4
    environment:
      POSTGRES_DB: geodata
      POSTGRES_USER: geouser
      POSTGRES_PASSWORD: ${POSTGRES_PASSWORD}
    volumes:
      - pgdata:/var/lib/postgresql/data
    ports:
      - "5432:5432"

  geoserver:
    image: kartoza/geoserver:2.24.2
    environment:
      GEOSERVER_DATA_DIR: /opt/geoserver/data_dir
      GEOSERVER_ADMIN_USER: admin
      GEOSERVER_ADMIN_PASSWORD: ${GEOSERVER_ADMIN_PASSWORD}
      INITIAL_MEMORY: 1G
      MAXIMUM_MEMORY: 4G
      STABLE_EXTENSIONS: "css-plugin,mbstyle-plugin"
      COMMUNITY_EXTENSIONS: "ogcapi-features-plugin"
    volumes:
      - geoserver_data:/opt/geoserver/data_dir
    ports:
      - "8080:8080"
    depends_on:
      - postgis

volumes:
  pgdata:
  geoserver_data:

REST API Configuration

# Create workspace
curl -u admin:geoserver -X POST \
  http://localhost:8080/geoserver/rest/workspaces \
  -H "Content-Type: application/json" \
  -d '{"workspace": {"name": "hamburg"}}'

# Create PostGIS data store
curl -u admin:geoserver -X POST \
  http://localhost:8080/geoserver/rest/workspaces/hamburg/datastores \
  -H "Content-Type: application/json" \
  -d '{
    "dataStore": {
      "name": "postgis_geodata",
      "connectionParameters": {
        "entry": [
          {"@key": "host", "$": "postgis"},
          {"@key": "port", "$": "5432"},
          {"@key": "database", "$": "geodata"},
          {"@key": "user", "$": "geouser"},
          {"@key": "passwd", "$": "secret"},
          {"@key": "dbtype", "$": "postgis"},
          {"@key": "schema", "$": "public"},
          {"@key": "Loose bbox", "$": "true"},
          {"@key": "Estimated extends", "$": "true"},
          {"@key": "Max open prepared statements", "$": "50"}
        ]
      }
    }
  }'

# Publish a layer from an existing PostGIS table
curl -u admin:geoserver -X POST \
  http://localhost:8080/geoserver/rest/workspaces/hamburg/datastores/postgis_geodata/featuretypes \
  -H "Content-Type: application/json" \
  -d '{
    "featureType": {
      "name": "facilities",
      "nativeName": "facilities",
      "title": "Healthcare Facilities",
      "srs": "EPSG:25832",
      "enabled": true
    }
  }'

SLD Styling

<?xml version="1.0" encoding="UTF-8"?>
<StyledLayerDescriptor version="1.0.0"
  xsi:schemaLocation="http://www.opengis.net/sld StyledLayerDescriptor.xsd"
  xmlns="http://www.opengis.net/sld"
  xmlns:ogc="http://www.opengis.net/ogc"
  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
  <NamedLayer>
    <Name>facilities_style</Name>
    <UserStyle>
      <Title>Healthcare Facilities</Title>
      <FeatureTypeStyle>
        <Rule>
          <Name>hospital</Name>
          <ogc:Filter>
            <ogc:PropertyIsEqualTo>
              <ogc:PropertyName>category</ogc:PropertyName>
              <ogc:Literal>hospital</ogc:Literal>
            </ogc:PropertyIsEqualTo>
          </ogc:Filter>
          <PointSymbolizer>
            <Graphic>
              <Mark>
                <WellKnownName>circle</WellKnownName>
                <Fill><CssParameter name="fill">#e74c3c</CssParameter></Fill>
                <Stroke>
                  <CssParameter name="stroke">#ffffff</CssParameter>
                  <CssParameter name="stroke-width">1</CssParameter>
                </Stroke>
              </Mark>
              <Size>12</Size>
            </Graphic>
          </PointSymbolizer>
        </Rule>
        <Rule>
          <Name>pharmacy</Name>
          <ogc:Filter>
            <ogc:PropertyIsEqualTo>
              <ogc:PropertyName>category</ogc:PropertyName>
              <ogc:Literal>pharmacy</ogc:Literal>
            </ogc:PropertyIsEqualTo>
          </ogc:Filter>
          <PointSymbolizer>
            <Graphic>
              <Mark>
                <WellKnownName>circle</WellKnownName>
                <Fill><CssParameter name="fill">#2ecc71</CssParameter></Fill>
                <Stroke>
                  <CssParameter name="stroke">#ffffff</CssParameter>
                  <CssParameter name="stroke-width">1</CssParameter>
                </Stroke>
              </Mark>
              <Size>8</Size>
            </Graphic>
          </PointSymbolizer>
        </Rule>
      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>

# Upload SLD style via REST API
curl -u admin:geoserver -X POST \
  http://localhost:8080/geoserver/rest/workspaces/hamburg/styles \
  -H "Content-Type: application/vnd.ogc.sld+xml" \
  -d @facilities_style.sld

# Assign style to layer
curl -u admin:geoserver -X PUT \
  http://localhost:8080/geoserver/rest/layers/hamburg:facilities \
  -H "Content-Type: application/json" \
  -d '{
    "layer": {
      "defaultStyle": {"name": "hamburg:facilities_style"}
    }
  }'

WMS Tile Serving

Raw WMS requests render images on the fly for every request. GeoWebCache (built into GeoServer) intercepts these requests, generates tiles once, and serves them from disk or memory on subsequent requests. This reduces server load by orders of magnitude and makes the difference between a sluggish map and a responsive one.

Tile Matrix Sets

Image Format Optimization

PNG supports transparency and is lossless -- ideal for overlays and vector-style maps. JPEG is lossy but produces smaller files -- suitable for satellite or aerial imagery where transparency is not needed. PNG8 (256 colours) is a useful middle ground for thematic maps with limited colour palettes, cutting file sizes by 60-70% compared to PNG24.

Cache Seeding

# Seed tiles via GeoWebCache REST API
curl -u admin:geoserver -X POST \
  http://localhost:8080/geoserver/gwc/rest/seed/hamburg:facilities.json \
  -H "Content-Type: application/json" \
  -d '{
    "seedRequest": {
      "name": "hamburg:facilities",
      "bounds": {
        "coords": {
          "double": [9.7, 53.4, 10.3, 53.7]
        }
      },
      "srs": {"number": 4326},
      "zoomStart": 8,
      "zoomStop": 16,
      "format": "image/png",
      "type": "seed",
      "threadCount": 4
    }
  }'

# Check seeding progress
curl -u admin:geoserver \
  http://localhost:8080/geoserver/gwc/rest/seed/hamburg:facilities.json

# Truncate cache for a layer (after data update)
curl -u admin:geoserver -X POST \
  http://localhost:8080/geoserver/gwc/rest/seed/hamburg:facilities.json \
  -H "Content-Type: application/json" \
  -d '{
    "seedRequest": {
      "name": "hamburg:facilities",
      "type": "truncate",
      "threadCount": 1
    }
  }'

Performance Optimization

Spatial queries are computationally expensive. A single ST_Intersection call on complex polygons can consume hundreds of milliseconds. Multiply that by thousands of features and you have a query that takes minutes. The optimizations below routinely reduce query times from minutes to sub-second.

Query Optimization with EXPLAIN ANALYZE

-- Always check the query plan for spatial queries
EXPLAIN (ANALYZE, BUFFERS, FORMAT TEXT)
SELECT d.name, COUNT(f.id)
FROM districts d
JOIN facilities f ON ST_Contains(d.geom, f.geom)
GROUP BY d.name;

-- Expected output should show:
-- Index Scan using idx_facilities_geom
-- Recheck Cond: (spatial filter)
-- Heap Blocks: exact=...

-- If you see Seq Scan on large tables, the index is not being used.
-- Common causes:
--   1. Missing spatial index
--   2. SRID mismatch between columns
--   3. Function wrapping that prevents index usage

Materialized Views for Complex Aggregations

-- Pre-compute facility counts per district
CREATE MATERIALIZED VIEW mv_district_facility_stats AS
SELECT
  d.id AS district_id,
  d.name,
  d.geom,
  COUNT(f.id) AS facility_count,
  ARRAY_AGG(DISTINCT f.category) AS categories,
  ST_Area(d.geom) / 1e6 AS area_km2,
  COUNT(f.id)::numeric / (ST_Area(d.geom) / 1e6) AS density_per_km2
FROM districts d
LEFT JOIN facilities f ON ST_Contains(d.geom, f.geom)
GROUP BY d.id, d.name, d.geom;

-- Add spatial index to the materialized view
CREATE INDEX idx_mv_district_stats_geom
  ON mv_district_facility_stats USING GIST (geom);

-- Refresh on schedule (e.g., nightly via pg_cron)
SELECT cron.schedule(
  'refresh_district_stats',
  '0 3 * * *',
  'REFRESH MATERIALIZED VIEW CONCURRENTLY mv_district_facility_stats'
);

Connection Pooling and Parallel Query

-- postgresql.conf tuning for spatial workloads
max_connections = 200
shared_buffers = 4GB
effective_cache_size = 12GB
work_mem = 256MB                 -- spatial joins need large work_mem
maintenance_work_mem = 1GB       -- for CREATE INDEX on large tables
max_parallel_workers_per_gather = 4
max_parallel_workers = 8
parallel_tuple_cost = 0.01       -- encourage parallel plans
random_page_cost = 1.1           -- SSD-optimized

-- PgBouncer config for GeoServer connection pooling
# pgbouncer.ini
[databases]
geodata = host=127.0.0.1 port=5432 dbname=geodata

[pgbouncer]
listen_port = 6432
listen_addr = 0.0.0.0
auth_type = md5
pool_mode = transaction
default_pool_size = 40
max_client_conn = 400
min_pool_size = 10

Integration with MapLibre GL

MapLibre GL JS renders interactive maps in the browser using WebGL. It can consume raster tiles from GeoServer WMS/WMTS endpoints and vector tiles from pg_tileserv, enabling both server-rendered and client-rendered map layers in the same application.

Consuming WMS Tiles

import maplibregl from "maplibre-gl";

const map = new maplibregl.Map({
  container: "map",
  style: {
    version: 8,
    sources: {
      osm: {
        type: "raster",
        tiles: ["https://tile.openstreetmap.org/{z}/{x}/{y}.png"],
        tileSize: 256,
        attribution: "OpenStreetMap contributors",
      },
    },
    layers: [{ id: "osm", type: "raster", source: "osm" }],
  },
  center: [9.9937, 53.5511],
  zoom: 11,
});

map.on("load", () => {
  // Add GeoServer WMS as raster source
  map.addSource("facilities-wms", {
    type: "raster",
    tiles: [
      "http://localhost:8080/geoserver/hamburg/wms?" +
        "service=WMS&version=1.1.1&request=GetMap" +
        "&layers=hamburg:facilities" +
        "&bbox={bbox-epsg-3857}" +
        "&width=256&height=256" +
        "&srs=EPSG:3857" +
        "&styles=facilities_style" +
        "&format=image/png&transparent=true",
    ],
    tileSize: 256,
  });

  map.addLayer({
    id: "facilities-layer",
    type: "raster",
    source: "facilities-wms",
    paint: { "raster-opacity": 0.85 },
  });
});

Vector Tiles with pg_tileserv

# Deploy pg_tileserv alongside PostGIS
# docker-compose addition:
  pg_tileserv:
    image: pramsey/pg_tileserv:latest
    environment:
      DATABASE_URL: postgresql://geouser:secret@postgis:5432/geodata
    ports:
      - "7800:7800"
    depends_on:
      - postgis

// Add vector tiles from pg_tileserv
map.addSource("districts-vector", {
  type: "vector",
  tiles: [
    "http://localhost:7800/public.districts/{z}/{x}/{y}.pbf",
  ],
  minzoom: 6,
  maxzoom: 16,
});

map.addLayer({
  id: "districts-fill",
  type: "fill",
  source: "districts-vector",
  "source-layer": "public.districts",
  paint: {
    "fill-color": [
      "interpolate",
      ["linear"],
      ["get", "population"],
      10000, "#f7fbff",
      50000, "#6baed6",
      100000, "#08306b",
    ],
    "fill-opacity": 0.6,
  },
});

map.addLayer({
  id: "districts-outline",
  type: "line",
  source: "districts-vector",
  "source-layer": "public.districts",
  paint: {
    "line-color": "#ffffff",
    "line-width": 1,
    "line-opacity": 0.4,
  },
});

Healthcare Use Case: Hamburg Isochrone Analysis

This section demonstrates a real-world application: analysing healthcare accessibility in Hamburg by generating walking isochrones around existing facilities, identifying coverage gaps, and optimising placement of new facilities using spatial analysis.

Data Preparation

-- Import Hamburg road network for routing
-- (Typically sourced from OpenStreetMap via osm2pgrouting)

-- Install pgRouting extension
CREATE EXTENSION IF NOT EXISTS pgrouting;

-- Load OSM data into routable network
-- osm2pgrouting command:
-- osm2pgrouting \
--   --f /data/osm/hamburg-latest.osm \
--   --conf /usr/share/osm2pgrouting/mapconfig_for_pedestrian.xml \
--   --dbname geodata \
--   --username geouser \
--   --clean

-- Add walking cost column (based on way length and speed)
ALTER TABLE road_network ADD COLUMN cost_walk DOUBLE PRECISION;
ALTER TABLE road_network ADD COLUMN reverse_cost_walk DOUBLE PRECISION;

UPDATE road_network
SET cost_walk = ST_Length(the_geom) / (5000.0 / 3600.0),  -- 5 km/h walking
    reverse_cost_walk = ST_Length(the_geom) / (5000.0 / 3600.0)
WHERE one_way != 1;  -- bidirectional

UPDATE road_network
SET reverse_cost_walk = 1e9
WHERE one_way = 1;  -- one-way streets

Isochrone Generation Pipeline

-- Generate isochrones for all healthcare facilities
CREATE TABLE healthcare_isochrones AS
WITH facility_vertices AS (
  -- Snap each facility to the nearest road network vertex
  SELECT
    f.id AS facility_id,
    f.name,
    (SELECT v.id
     FROM road_network_vertices_pgr v
     ORDER BY v.the_geom <-> f.geom
     LIMIT 1) AS vertex_id
  FROM facilities f
  WHERE f.category = 'healthcare'
),
isochrone_nodes AS (
  SELECT
    fv.facility_id,
    fv.name,
    dd.agg_cost,
    v.the_geom
  FROM facility_vertices fv,
  LATERAL pgr_drivingDistance(
    'SELECT id, source, target, cost_walk AS cost,
            reverse_cost_walk AS reverse_cost
     FROM road_network',
    fv.vertex_id,
    2700,  -- max 45 minutes
    directed := false
  ) dd
  JOIN road_network_vertices_pgr v ON dd.node = v.id
)
SELECT
  facility_id,
  name,
  CASE
    WHEN agg_cost <= 900 THEN 15
    WHEN agg_cost <= 1800 THEN 30
    ELSE 45
  END AS minutes,
  ST_Transform(
    ST_ConcaveHull(ST_Collect(the_geom), 0.7),
    25832
  ) AS geom
FROM isochrone_nodes
GROUP BY facility_id, name,
  CASE
    WHEN agg_cost <= 900 THEN 15
    WHEN agg_cost <= 1800 THEN 30
    ELSE 45
  END;

-- Add spatial index
CREATE INDEX idx_isochrones_geom
  ON healthcare_isochrones USING GIST (geom);

Coverage Gap Identification

-- Find residential areas outside 30-minute walking coverage
WITH coverage_30min AS (
  SELECT ST_Union(geom) AS covered
  FROM healthcare_isochrones
  WHERE minutes <= 30
),
residential_areas AS (
  SELECT ST_Union(geom) AS residential
  FROM buildings
  WHERE building_type = 'residential'
)
SELECT
  ST_Intersection(ra.residential, ST_Difference(
    (SELECT ST_Union(geom) FROM districts),
    c30.covered
  )) AS underserved_residential,
  ST_Area(ST_Intersection(ra.residential, ST_Difference(
    (SELECT ST_Union(geom) FROM districts),
    c30.covered
  ))) / 1e6 AS underserved_km2
FROM coverage_30min c30, residential_areas ra;

-- Rank potential locations for new facilities
-- by population density in underserved areas
WITH gap_zones AS (
  SELECT ST_Difference(
    (SELECT ST_Union(geom) FROM districts),
    (SELECT ST_Union(geom) FROM healthcare_isochrones WHERE minutes <= 30)
  ) AS gap_geom
),
candidate_grid AS (
  SELECT (ST_SquareGrid(500, gap_geom)).*  -- 500m grid cells
  FROM gap_zones
)
SELECT
  cg.geom AS grid_cell,
  COUNT(b.id) AS building_count,
  SUM(b.estimated_residents) AS population_estimate
FROM candidate_grid cg
JOIN buildings b ON ST_Intersects(cg.geom, b.geom)
WHERE b.building_type = 'residential'
GROUP BY cg.geom
ORDER BY population_estimate DESC
LIMIT 20;

Production Operations

Running PostGIS and GeoServer in production requires attention to backup strategies, monitoring, automated deployments, and schema migration. The patterns below are drawn from operating spatial infrastructure serving millions of tile requests per day.

Backup Strategies for Spatial Data

#!/bin/bash
# backup_geodata.sh -- PostGIS backup with spatial data preservation

set -euo pipefail

BACKUP_DIR="/backups/postgis/$(date +%Y-%m-%d)"
DB_NAME="geodata"
DB_USER="geouser"

mkdir -p "$BACKUP_DIR"

# Full custom-format dump (preserves PostGIS types and indexes)
pg_dump -U "$DB_USER" -Fc -Z 6 \
  --no-owner --no-acl \
  -f "$BACKUP_DIR/${DB_NAME}_full.dump" \
  "$DB_NAME"

# Schema-only dump for version control
pg_dump -U "$DB_USER" --schema-only \
  -f "$BACKUP_DIR/${DB_NAME}_schema.sql" \
  "$DB_NAME"

# Export critical layers as GeoPackage (portable format)
ogr2ogr -f "GPKG" \
  "$BACKUP_DIR/facilities.gpkg" \
  PG:"host=localhost dbname=$DB_NAME user=$DB_USER" \
  facilities districts healthcare_isochrones

# Rotate: keep last 30 days
find /backups/postgis/ -maxdepth 1 -mtime +30 -type d -exec rm -rf {} +

echo "Backup complete: $BACKUP_DIR"

GeoServer Health Monitoring

# Prometheus metrics from GeoServer monitoring module
# Enable in geoserver/WEB-INF/web.xml:
#   <context-param>
#     <param-name>GEOSERVER_CSRF_DISABLED</param-name>
#     <param-value>true</param-value>
#   </context-param>

# Health check script
#!/bin/bash
set -euo pipefail

GEOSERVER_URL="http://localhost:8080/geoserver"

# Check GeoServer is responding
HTTP_CODE=$(curl -s -o /dev/null -w "%{http_code}" \
  "$GEOSERVER_URL/web/")

if [ "$HTTP_CODE" != "200" ]; then
  echo "CRITICAL: GeoServer returned HTTP $HTTP_CODE"
  exit 2
fi

# Check WMS GetCapabilities
WMS_CODE=$(curl -s -o /dev/null -w "%{http_code}" \
  "$GEOSERVER_URL/hamburg/wms?service=WMS&request=GetCapabilities")

if [ "$WMS_CODE" != "200" ]; then
  echo "WARNING: WMS GetCapabilities returned HTTP $WMS_CODE"
  exit 1
fi

# Check tile cache disk usage
CACHE_SIZE=$(du -sb /opt/geoserver/data_dir/gwc/ | cut -f1)
CACHE_GB=$((CACHE_SIZE / 1073741824))

if [ "$CACHE_GB" -gt 50 ]; then
  echo "WARNING: GeoWebCache using ${CACHE_GB}GB"
  exit 1
fi

echo "OK: GeoServer healthy, cache ${CACHE_GB}GB"
exit 0

CI/CD for Spatial Schemas

# .github/workflows/spatial-deploy.yml
name: Deploy Spatial Schema

on:
  push:
    branches: [main]
    paths:
      - "sql/migrations/**"
      - "styles/**"

jobs:
  deploy:
    runs-on: ubuntu-latest
    steps:
      - uses: actions/checkout@v4

      - name: Run PostGIS migrations
        env:
          PGHOST: ${{ secrets.POSTGIS_HOST }}
          PGDATABASE: geodata
          PGUSER: ${{ secrets.POSTGIS_USER }}
          PGPASSWORD: ${{ secrets.POSTGIS_PASSWORD }}
        run: |
          for f in sql/migrations/*.sql; do
            echo "Applying: $f"
            psql -v ON_ERROR_STOP=1 -f "$f"
          done

      - name: Deploy SLD styles to GeoServer
        run: |
          for style in styles/*.sld; do
            STYLE_NAME=$(basename "$style" .sld)
            curl -u admin:${{ secrets.GEOSERVER_PASSWORD }} \
              -X PUT \
              "${{ secrets.GEOSERVER_URL }}/rest/workspaces/hamburg/styles/$STYLE_NAME" \
              -H "Content-Type: application/vnd.ogc.sld+xml" \
              -d @"$style"
          done

      - name: Truncate affected tile caches
        run: |
          curl -u admin:${{ secrets.GEOSERVER_PASSWORD }} \
            -X POST \
            "${{ secrets.GEOSERVER_URL }}/gwc/rest/masstruncate" \
            -H "Content-Type: application/json" \
            -d '{"truncateLayer": {"layerName": "hamburg:facilities"}}'

      - name: Refresh materialized views
        env:
          PGHOST: ${{ secrets.POSTGIS_HOST }}
          PGDATABASE: geodata
          PGUSER: ${{ secrets.POSTGIS_USER }}
          PGPASSWORD: ${{ secrets.POSTGIS_PASSWORD }}
        run: |
          psql -c "REFRESH MATERIALIZED VIEW CONCURRENTLY mv_district_facility_stats;"

Automated Style Deployment

#!/usr/bin/env python3
"""deploy_styles.py -- Sync SLD styles to GeoServer via REST API."""

import os
import sys
from pathlib import Path
import requests

GEOSERVER_URL = os.environ["GEOSERVER_URL"]
GEOSERVER_USER = os.environ.get("GEOSERVER_USER", "admin")
GEOSERVER_PASS = os.environ["GEOSERVER_PASSWORD"]
WORKSPACE = os.environ.get("GEOSERVER_WORKSPACE", "hamburg")

STYLES_DIR = Path("styles")
session = requests.Session()
session.auth = (GEOSERVER_USER, GEOSERVER_PASS)


def deploy_style(style_path: Path) -> bool:
    """Upload or update a single SLD style."""
    style_name = style_path.stem
    sld_content = style_path.read_text(encoding="utf-8")

    url = f"{GEOSERVER_URL}/rest/workspaces/{WORKSPACE}/styles/{style_name}"

    # Check if style exists
    resp = session.get(url)
    if resp.status_code == 200:
        # Update existing
        resp = session.put(
            url,
            data=sld_content,
            headers={"Content-Type": "application/vnd.ogc.sld+xml"},
        )
    else:
        # Create new
        resp = session.post(
            f"{GEOSERVER_URL}/rest/workspaces/{WORKSPACE}/styles",
            data=sld_content,
            headers={"Content-Type": "application/vnd.ogc.sld+xml"},
            params={"name": style_name},
        )

    if resp.ok:
        print(f"  Deployed: {style_name}")
        return True
    else:
        print(f"  FAILED: {style_name} -- {resp.status_code} {resp.text}")
        return False


def main():
    sld_files = sorted(STYLES_DIR.glob("*.sld"))
    print(f"Found {len(sld_files)} styles to deploy.")

    results = [deploy_style(f) for f in sld_files]
    failed = results.count(False)

    if failed:
        print(f"\n{failed} style(s) failed to deploy.")
        sys.exit(1)

    print("\nAll styles deployed successfully.")


if __name__ == "__main__":
    main()

Conclusion

PostGIS and GeoServer remain the most capable open-source combination for spatial data infrastructure. PostGIS gives you a SQL-native spatial engine that integrates with everything PostgreSQL offers -- ACID transactions, MVCC, logical replication, and a mature ecosystem of extensions like pgRouting and pg_tileserv. GeoServer adds standardised OGC service endpoints and a tile cache that can serve millions of requests per day from commodity hardware.

The key architectural decisions that determine success are: choosing the right SRID and geometry type upfront, building spatial indexes before loading large datasets, using materialized views for complex aggregations that feed map layers, and configuring GeoWebCache to absorb read traffic. With these foundations in place, you can build spatial applications that scale from prototype to production without rewriting the data layer.