Serving GeoJSON, vector tiles and OGC APIs

geojson_session.py

"""GeoJSON from SQL: a point as stored and as RFC 7946 wants it, a row as a Feature, and whole collections."""
from pg_session import Session

with Session() as pg:
    pg.run("python3 make_town.py")
    pg.run("createdb town")
    pg.run("psql -X -q -d town -f load_town.sql")
    pg.psql("geojson.sql", db="town")

$ python3 make_town.py
wrote wards.csv (12 rows), households.csv (5,000 rows) and clinics.csv (8 rows)

$ createdb town

$ psql -X -q -d town -f load_town.sql
 wards | households | clinics 
-------+------------+---------
    12 |       5000 |       8
(1 row)

$ psql -X -a -d town -f geojson.sql
-- Clinic C2 as stored: metres of UTM zone 43N, and a crs member that RFC 7946 no longer has.
SELECT ST_AsGeoJSON(geom) AS as_stored FROM clinic WHERE clinic = 'C2';
                                                as_stored                                                 
----------------------------------------------------------------------------------------------------------
 {"type":"Point","crs":{"type":"name","properties":{"name":"EPSG:32643"}},"coordinates":[504000,2005000]}
(1 row)

-- Transformed to longitude and latitude, rounded to 6 decimal places.
SELECT ST_AsGeoJSON(ST_Transform(geom, 4326), 6) AS for_the_web FROM clinic WHERE clinic = 'C2';
                     for_the_web                      
------------------------------------------------------
 {"type":"Point","coordinates":[75.037813,18.133898]}
(1 row)

-- A whole row as a Feature: every other column becomes a property.
SELECT ST_AsGeoJSON(t.*, maxdecimaldigits => 6) AS feature
FROM (SELECT clinic, ST_Transform(geom, 4326) AS geom FROM clinic WHERE clinic = 'C2') AS t;
                                                        feature                                                        
-----------------------------------------------------------------------------------------------------------------------
 {"type": "Feature", "geometry": {"type":"Point","coordinates":[75.037813,18.133898]}, "properties": {"clinic": "C2"}}
(1 row)

-- The 8 clinics as one FeatureCollection, assembled in the query.
SELECT json_build_object('type', 'FeatureCollection', 'features', json_agg(ST_AsGeoJSON(t.*, maxdecimaldigits => 6)::json))
         ->'features'->0 AS first_of_8,
       length(json_build_object('type', 'FeatureCollection', 'features',
                                json_agg(ST_AsGeoJSON(t.*, maxdecimaldigits => 6)::json))::text) AS bytes
FROM (SELECT clinic, ST_Transform(geom, 4326) AS geom FROM clinic ORDER BY clinic) AS t;
                                                      first_of_8                                                      | bytes 
----------------------------------------------------------------------------------------------------------------------+-------
 {"type": "Feature", "geometry": {"type":"Point","coordinates":[75.01418,18.135709]}, "properties": {"clinic": "C1"}} |   995
(1 row)

-- All 5,000 households: the size of the collection at 15, 9 and 6 decimal places.
SELECT d AS decimals, length(json_build_object('type', 'FeatureCollection', 'features',
         json_agg(ST_AsGeoJSON(t.*, maxdecimaldigits => d)::json))::text) AS bytes
FROM (SELECT household, ST_Transform(geom, 4326) AS geom FROM household) AS t, (VALUES (15), (9), (6)) AS v (d)
GROUP BY d ORDER BY d DESC;
 decimals | bytes  
----------+--------
       15 | 707344
        9 | 653954
        6 | 623883
(3 rows)

geojson.sql

-- Clinic C2 as stored: metres of UTM zone 43N, and a crs member that RFC 7946 no longer has.
SELECT ST_AsGeoJSON(geom) AS as_stored FROM clinic WHERE clinic = 'C2';
-- Transformed to longitude and latitude, rounded to 6 decimal places.
SELECT ST_AsGeoJSON(ST_Transform(geom, 4326), 6) AS for_the_web FROM clinic WHERE clinic = 'C2';
-- A whole row as a Feature: every other column becomes a property.
SELECT ST_AsGeoJSON(t.*, maxdecimaldigits => 6) AS feature
FROM (SELECT clinic, ST_Transform(geom, 4326) AS geom FROM clinic WHERE clinic = 'C2') AS t;

-- The 8 clinics as one FeatureCollection, assembled in the query.
SELECT json_build_object('type', 'FeatureCollection', 'features', json_agg(ST_AsGeoJSON(t.*, maxdecimaldigits => 6)::json))
         ->'features'->0 AS first_of_8,
       length(json_build_object('type', 'FeatureCollection', 'features',
                                json_agg(ST_AsGeoJSON(t.*, maxdecimaldigits => 6)::json))::text) AS bytes
FROM (SELECT clinic, ST_Transform(geom, 4326) AS geom FROM clinic ORDER BY clinic) AS t;

-- All 5,000 households: the size of the collection at 15, 9 and 6 decimal places.
SELECT d AS decimals, length(json_build_object('type', 'FeatureCollection', 'features',
         json_agg(ST_AsGeoJSON(t.*, maxdecimaldigits => d)::json))::text) AS bytes
FROM (SELECT household, ST_Transform(geom, 4326) AS geom FROM household) AS t, (VALUES (15), (9), (6)) AS v (d)
GROUP BY d ORDER BY d DESC;

load_town.sql

-- The town in three typed tables, from the CSV files make_town.py writes (metres of UTM zone 43N).
CREATE EXTENSION postgis;
CREATE TABLE ward (ward text PRIMARY KEY, geom geometry(Polygon, 32643) NOT NULL);
CREATE TABLE household (household text PRIMARY KEY, geom geometry(Point, 32643) NOT NULL);
CREATE TABLE clinic (clinic text PRIMARY KEY, geom geometry(Point, 32643) NOT NULL);

CREATE TEMP TABLE ward_csv (ward text, wkt text);
\copy ward_csv FROM 'wards.csv' WITH (FORMAT csv, HEADER true)
INSERT INTO ward SELECT ward, ST_GeomFromText(wkt, 32643) FROM ward_csv;

CREATE TEMP TABLE point_csv (id text, x float8, y float8);
\copy point_csv FROM 'households.csv' WITH (FORMAT csv, HEADER true)
INSERT INTO household SELECT id, ST_Point(x, y, 32643) FROM point_csv;
TRUNCATE point_csv;
\copy point_csv FROM 'clinics.csv' WITH (FORMAT csv, HEADER true)
INSERT INTO clinic SELECT id, ST_Point(x, y, 32643) FROM point_csv;

CREATE INDEX ON ward USING gist (geom);
CREATE INDEX ON household USING gist (geom);
CREATE INDEX ON clinic USING gist (geom);
ANALYZE ward, household, clinic;
SELECT (SELECT count(*) FROM ward) AS wards, (SELECT count(*) FROM household) AS households,
       (SELECT count(*) FROM clinic) AS clinics;

make_town.py

"""Writes the track's town as three CSV files in metres of UTM zone 43N (EPSG:32643): the 12 wards as WKT polygons,
the 5,000 households and the 8 clinics as x and y. Python's standard library only."""
import csv

from town import clinics, households, ward_bounds

with open("wards.csv", "w", newline="") as f:
    out = csv.writer(f)
    out.writerow(["ward", "wkt"])
    for ward, (x0, y0, x1, y1) in ward_bounds().items():
        out.writerow([ward, f"POLYGON (({x0} {y0}, {x1} {y0}, {x1} {y1}, {x0} {y1}, {x0} {y0}))"])

for name, rows in (("households", households()), ("clinics", clinics())):
    with open(f"{name}.csv", "w", newline="") as f:
        out = csv.writer(f)
        out.writerow([name[:-1], "x", "y"])
        out.writerows(rows)

print("wrote wards.csv (12 rows), households.csv (5,000 rows) and clinics.csv (8 rows)")

pg_session.py

"""Runs one session of this lesson on a throwaway PostgreSQL 18.6 server with PostGIS 3.6.4.

The session makes a new cluster in a folder of its own, starts it on 127.0.0.1 only (a port the operating system
reports as free, no Unix socket, a password made up for the run), runs psql scripts and the command-line tools of
PostGIS and GDAL against it, and prints each command with what it printed. Then it stops the server and deletes the
cluster, also when something fails. The tools read the connection from the standard libpq variables (PGHOST, PGPORT,
PGUSER, PGPASSWORD), so no command shows the port or the password.

The tools are looked up in the folder that PG_BIN names, then in the bin folder of a folder named postgis-toolchain in
the session's folder or in any folder above it, then on the PATH. To repeat a session on your own computer, install
PostgreSQL 18 with PostGIS 3.6 and GDAL 3.13, put their bin folder on the PATH and type the commands after "$ ".

The output is what the programs printed, with the session's folder shown as "." so that it reads the same on every
computer. A command that should succeed and fails, or a script that prints another number of ERROR lines than the
session expects, ends the session with status 1, so a run that went wrong is never kept as if it had worked.
"""

import os
import secrets
import shutil
import socket
import subprocess
import sys
import tempfile
from pathlib import Path


def _tool_folder():
    """The folder that holds postgres, psql and the other tools."""
    named = os.environ.get("PG_BIN")
    if named:
        return Path(named)
    here = Path.cwd().resolve()
    for folder in (here, *here.parents):
        candidate = folder / "postgis-toolchain" / "bin"
        if (candidate / "postgres").exists():
            return candidate
    found = shutil.which("postgres")
    if found:
        return Path(found).resolve().parent
    sys.exit("PostgreSQL with PostGIS was not found: set PG_BIN or put its bin folder on the PATH")


BIN = _tool_folder()
# A pipeline (shp2pgsql … | psql …) fails when any of its commands fails: bash's pipefail.
SHELL = shutil.which("bash") or "/bin/sh"

# Settings that make every run print the same: no parallel plans, no just-in-time compilation, no background vacuum.
SETTINGS = {
    "jit": "off",
    "max_parallel_workers_per_gather": "0",
    "autovacuum": "off",
    "timezone": "'UTC'",
    "lc_messages": "'C'",
}


def _free_port():
    """A TCP port on 127.0.0.1 that nothing listens on now."""
    with socket.socket() as probe:
        probe.bind(("127.0.0.1", 0))
        return probe.getsockname()[1]


class Session:
    """A cluster that exists for one with block: made and started on entry, stopped and deleted on exit."""

    def __init__(self, settings=None):
        self.settings = {**SETTINGS, **(settings or {})}
        self.failed = False

    def __enter__(self):
        self.folder = Path.cwd()
        self.base = Path(tempfile.mkdtemp(prefix="cluster-", dir=self.folder))
        self.data = self.base / "data"
        password = secrets.token_urlsafe(24)
        (self.base / "pwfile").write_text(password + "\n")
        self.env = {**os.environ, "PATH": f"{BIN}{os.pathsep}{os.environ.get('PATH', '')}", "PGHOST": "127.0.0.1",
                    "PGPORT": str(_free_port()), "PGUSER": "postgres", "PGPASSWORD": password,
                    "PGDATABASE": "postgres", "PGCONNECT_TIMEOUT": "10"}
        for key in ("PGSERVICE", "PGSERVICEFILE", "PGOPTIONS", "PSQLRC", "PGSSLMODE"):
            self.env.pop(key, None)
        self._quiet(["initdb", "-D", str(self.data), "-U", "postgres", "--auth=scram-sha-256",
                     f"--pwfile={self.base / 'pwfile'}", "--no-locale", "--encoding=UTF8", "--no-instructions"])
        with open(self.data / "postgresql.conf", "a") as conf:
            conf.write(f"listen_addresses = '127.0.0.1'\nport = {self.env['PGPORT']}\nunix_socket_directories = ''\n")
            for key, value in self.settings.items():
                conf.write(f"{key} = {value}\n")
        self._quiet(["pg_ctl", "-D", str(self.data), "-l", str(self.base / "server.log"), "-w", "-t", "60", "start"])
        return self

    def __exit__(self, *exc):
        try:
            self._quiet(["pg_ctl", "-D", str(self.data), "-m", "fast", "-w", "-t", "60", "stop"])
        finally:
            shutil.rmtree(self.base, ignore_errors=True)
        if self.failed and exc[0] is None:
            sys.exit(1)

    def _quiet(self, args):
        """Runs one of the server's own programs without printing it; stops the session if it fails."""
        done = subprocess.run([str(BIN / args[0]), *args[1:]], env=self.env, capture_output=True, text=True)
        if done.returncode != 0:
            sys.exit(f"{args[0]} failed: {done.stdout}{done.stderr}")

    def _clean(self, text):
        """The output as it reads on any computer: the session's folder as "." and no folder of the tools."""
        for folder in {str(self.folder), str(self.folder.resolve())}:
            text = text.replace(folder + "/", "").replace(folder, ".")
        return text.replace(str(BIN) + "/", "")

    def run(self, command, status=0, note=None):
        """Runs a command line as a terminal would, prints it and its output, and checks how it ended."""
        if note:
            print(f"# {note}")
        print(f"$ {command}")
        done = subprocess.run([SHELL, "-c", f"set -o pipefail; {command}"], cwd=self.folder, env=self.env,
                              stdout=subprocess.PIPE, stderr=subprocess.STDOUT, text=True)
        out = self._clean(done.stdout)
        if out.strip():
            print(out.rstrip("\n"))
        if done.returncode != status:
            print(f"[exit status {done.returncode}; the session expected {status}]")
            self.failed = True
        print()
        return out

    def psql(self, script, db="postgres", errors=0, note=None):
        """Runs a psql script with every line echoed (-a) and checks how many ERROR lines it printed."""
        out = self.run(f"psql -X -a -d {db} -f {script}", note=note)
        found = sum(1 for line in out.splitlines() if "ERROR:" in line)
        if found != errors:
            print(f"[{script}: {found} errors; the session expected {errors}]")
            self.failed = True
        return out

tile_session.py

"""One vector tile from SQL: the bounds of tile 12/2901/1838, its households encoded with ST_AsMVT, and one
household's place in tile units."""
from pg_session import Session

with Session() as pg:
    pg.run("python3 make_town.py")
    pg.run("createdb town")
    pg.run("psql -X -q -d town -f load_town.sql")
    pg.psql("tile.sql", db="town")

$ python3 make_town.py
wrote wards.csv (12 rows), households.csv (5,000 rows) and clinics.csv (8 rows)

$ createdb town

$ psql -X -q -d town -f load_town.sql
 wards | households | clinics 
-------+------------+---------
    12 |       5000 |       8
(1 row)

$ psql -X -a -d town -f tile.sql
-- The zoom-12 tile that holds clinic C2: column 2901, row 1838. Its bounds, in Web Mercator metres.
SELECT ST_AsText(ST_TileEnvelope(12, 2901, 1838), 0) AS bounds_3857;
                                        bounds_3857                                         
--------------------------------------------------------------------------------------------
 POLYGON((8345700 2044843,8345700 2054627,8355484 2054627,8355484 2044843,8345700 2044843))
(1 row)

-- The households of that tile in tile coordinates, encoded as one layer of a vector tile.
-- The tile's box is transformed into the data's SRID once, so the index on household.geom can answer.
WITH bounds AS (SELECT ST_TileEnvelope(12, 2901, 1838) AS b),
rows AS (
  SELECT h.household, ST_AsMVTGeom(ST_Transform(h.geom, 3857), bounds.b) AS geom
  FROM household h, bounds
  WHERE ST_Intersects(h.geom, ST_Transform(bounds.b, 32643)))
SELECT count(*) AS features, length(ST_AsMVT(rows.*, 'households')) AS tile_bytes FROM rows;
 features | tile_bytes 
----------+------------
     3288 |      82104
(1 row)

-- One household's position in tile units: 0 to 4096 across the tile, measured from its top-left corner.
SELECT h.household, ST_AsText(ST_AsMVTGeom(ST_Transform(h.geom, 3857), ST_TileEnvelope(12, 2901, 1838))) AS in_tile
FROM household h WHERE h.household = 'H0001';
 household |     in_tile     
-----------+-----------------
 H0001     | POINT(2280 966)
(1 row)

tile.sql

-- The zoom-12 tile that holds clinic C2: column 2901, row 1838. Its bounds, in Web Mercator metres.
SELECT ST_AsText(ST_TileEnvelope(12, 2901, 1838), 0) AS bounds_3857;

-- The households of that tile in tile coordinates, encoded as one layer of a vector tile.
-- The tile's box is transformed into the data's SRID once, so the index on household.geom can answer.
WITH bounds AS (SELECT ST_TileEnvelope(12, 2901, 1838) AS b),
rows AS (
  SELECT h.household, ST_AsMVTGeom(ST_Transform(h.geom, 3857), bounds.b) AS geom
  FROM household h, bounds
  WHERE ST_Intersects(h.geom, ST_Transform(bounds.b, 32643)))
SELECT count(*) AS features, length(ST_AsMVT(rows.*, 'households')) AS tile_bytes FROM rows;

-- One household's position in tile units: 0 to 4096 across the tile, measured from its top-left corner.
SELECT h.household, ST_AsText(ST_AsMVTGeom(ST_Transform(h.geom, 3857), ST_TileEnvelope(12, 2901, 1838))) AS in_tile
FROM household h WHERE h.household = 'H0001';