GIS & Geospatial Data (QGIS, GeoPandas, PostGIS) Module 8 – Spatial SQL with PostGIS

Geometry processing in PostGIS

Repair, buffer and union geometries in SQL, simplify shared borders without gaps, cluster points with DBSCAN and k-means, and count them on grids.

Advanced
15 minutes
Examples run with Python 3.14.8, Pyodide 314.0.7, PostgreSQL 18.6 and postgis 3.6.4
By MySmartCoPilot

Free: Repair first: name the fault, then keep what was drawn · Buffers and unions work in the units of the SRID. With a pass: 2 in-depth parts and 3 interview answers, about 9 minutes.

What you will learn

Repair invalid geometries with ST_MakeValid and check that the repair kept the intended shape
Build buffers and unions in a projected CRS, and count what they cover with ST_DWithin instead
Simplify a coverage of districts that share borders without opening gaps or overlaps
Group points into clusters with ST_ClusterDBSCAN and ST_ClusterKMeans, and count them on grids and geohash cells
Decide which processing belongs in the database and which in Python or QGIS

Before you start

Spatial relationships and joins in SQL

What you’ll learn · Before you start

What you will learn

Repair invalid geometries with ST_MakeValid and check that the repair kept the intended shape
Build buffers and unions in a projected CRS, and count what they cover with ST_DWithin instead
Simplify a coverage of districts that share borders without opening gaps or overlaps
Group points into clusters with ST_ClusterDBSCAN and ST_ClusterKMeans, and count them on grids and geohash cells
Decide which processing belongs in the database and which in Python or QGIS

Before you start

Spatial relationships and joins in SQL

On this page

A bow-tie polygon, one ring that crosses itself, reports an area of 0, because its two triangles cancel each other out. ST_Buffer(geom, 0), the repair many scripts still use, returns one triangle and silently loses the other; ST_MakeValid keeps both.

Here you repair, buffer and merge geometry where it is stored, in PostgreSQL 18.6 with PostGIS 3.6.4, and check each result.

Repair first: name the fault, then keep what was drawn

Every overlay, area and simplification in this lesson assumes valid geometry: A polygon or multipolygon that follows the simple-features rules: simple rings that do not cross, touch only at points and keep the holes inside the shell, and an interior in one piece. Areas and overlays are reliable only on valid geometry. More in “Geometry types, WKT, WKB and validity”, so the first step of any processing is to find the invalid rows and say why they are invalid. ST_IsValidReason returns the reason and the place, and ST_MakeValid builds a valid geometry “without losing any of the input vertices” (ST_MakeValid). Here is what the two repairs make of one bow-tie:

One bow-tie ring and two repairs

Text description of the diagram

Left, as written: a ring through the points 0 0, 2 2, 2 0 and 0 2, numbered 1 to 4 in the order the ring visits them, back to 0 0. Its two diagonals cross at 1 1, marked with a ring. The two triangles it encloses have areas of opposite sign, so ST_Area returns 0.

Middle, ST_MakeValid: a MultiPolygon of the two triangles, one on the left from 0 0 to 0 2 with its tip at 1 1, one on the right from 2 0 to 2 2 with its tip at 1 1. Area 2.

Right, ST_Buffer(g, 0): only the right-hand triangle is kept, drawn filled; the left one, drawn dashed, is lost. Area 1.

Step through how each repair rebuilds the ring:

Interactive · stage by stage The bow-tie through a validity check and two repairs Step through the checks and the repair, then follow the branch that buffers by zero.

The drawing shows the whole scene; its 7 steps are listed below.

Every step as text

The ring. The ring visits 0 0, 2 2, 2 0, 0 2 and closes at 0 0, so its two diagonals cross. ST_Area returns 0: the two triangles it encloses have areas of opposite sign.
The check. ST_IsValid is false and ST_IsValidReason says Self-intersection[1 1]: the fault and the place, which a report of invalid rows should keep.
Lines first. The default repair, linework, does not trust the ring. It takes its four edges as plain lines, whatever order they were drawn in.
Noding. The lines are split where they cross: both diagonals at 1 1. Every vertex of the input is kept, and 1 1 becomes a vertex of its own.
Two triangles. The noded lines close two loops, a triangle on each side of 1 1, and each becomes a polygon.
The result. ST_MakeValid returns MULTIPOLYGON(((0 2,1 1,0 0,0 2)),((2 0,1 1,2 2,2 0))), area 2: both halves of what was drawn.
The zero buffer. ST_Buffer(g, 0) computes a buffer of width zero around the ring as drawn and returns POLYGON((1 1,2 2,2 0,1 1)), area 1. The left triangle is gone, and nothing says so.

A bow-tie, two repairs, and a ward digitised with two corners swapped

repair_session.py

"""What ST_IsValidReason reports for a ring that crosses itself, and what two repairs make of it."""
from pg_session import Session

with Session() as pg:
    pg.run("createdb town")
    pg.run('psql -X -q -d town -c "CREATE EXTENSION postgis"')
    pg.psql("repair.sql", db="town")

Output

$ createdb town

$ psql -X -q -d town -c "CREATE EXTENSION postgis"

$ psql -X -a -d town -f repair.sql
-- A bow-tie: one ring that crosses itself at (1 1). Its two triangles have areas of opposite sign.
SELECT ST_IsValid(g) AS valid, ST_IsValidReason(g) AS reason, ST_Area(g) AS area
FROM (VALUES ('POLYGON((0 0, 2 2, 2 0, 0 2, 0 0))'::geometry)) AS t (g);
psql:repair.sql:3: NOTICE:  Self-intersection at or near point 1 1
 valid |         reason         | area 
-------+------------------------+------
 f     | Self-intersection[1 1] |    0
(1 row)

-- ST_MakeValid keeps both triangles; the old habit of buffering by zero keeps one.
SELECT 'ST_MakeValid' AS repair, ST_AsText(ST_MakeValid(g)) AS result, ST_Area(ST_MakeValid(g)) AS area
FROM (VALUES ('POLYGON((0 0, 2 2, 2 0, 0 2, 0 0))'::geometry)) AS t (g)
UNION ALL
SELECT 'ST_Buffer(g, 0)', ST_AsText(ST_Buffer(g, 0)), ST_Area(ST_Buffer(g, 0))
FROM (VALUES ('POLYGON((0 0, 2 2, 2 0, 0 2, 0 0))'::geometry)) AS t (g);
     repair      |                        result                         | area 
-----------------+-------------------------------------------------------+------
 ST_MakeValid    | MULTIPOLYGON(((0 2,1 1,0 0,0 2)),((2 0,1 1,2 2,2 0))) |    2
 ST_Buffer(g, 0) | POLYGON((1 1,2 2,2 0,1 1))                            |    1
(2 rows)

-- Ward W02 digitised with two corners swapped: repaired, it is two triangles, not the 4 km² square.
SELECT ST_IsValidReason(g) AS reason, ST_Area(g) AS area_m2, ST_Area(ST_MakeValid(g)) AS repaired_m2,
       ST_NumGeometries(ST_MakeValid(g)) AS parts
FROM (VALUES (ST_GeomFromText('POLYGON((502000 2004000, 504000 2004000, 502000 2006000, 504000 2006000,
                                        502000 2004000))', 32643))) AS t (g);
              reason               | area_m2 | repaired_m2 | parts 
-----------------------------------+---------+-------------+-------
 Self-intersection[503000 2005000] |       0 |     2000000 |     2
(1 row)

Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 repair_session.py

repair.sql

-- A bow-tie: one ring that crosses itself at (1 1). Its two triangles have areas of opposite sign.
SELECT ST_IsValid(g) AS valid, ST_IsValidReason(g) AS reason, ST_Area(g) AS area
FROM (VALUES ('POLYGON((0 0, 2 2, 2 0, 0 2, 0 0))'::geometry)) AS t (g);

-- ST_MakeValid keeps both triangles; the old habit of buffering by zero keeps one.
SELECT 'ST_MakeValid' AS repair, ST_AsText(ST_MakeValid(g)) AS result, ST_Area(ST_MakeValid(g)) AS area
FROM (VALUES ('POLYGON((0 0, 2 2, 2 0, 0 2, 0 0))'::geometry)) AS t (g)
UNION ALL
SELECT 'ST_Buffer(g, 0)', ST_AsText(ST_Buffer(g, 0)), ST_Area(ST_Buffer(g, 0))
FROM (VALUES ('POLYGON((0 0, 2 2, 2 0, 0 2, 0 0))'::geometry)) AS t (g);

-- Ward W02 digitised with two corners swapped: repaired, it is two triangles, not the 4 km² square.
SELECT ST_IsValidReason(g) AS reason, ST_Area(g) AS area_m2, ST_Area(ST_MakeValid(g)) AS repaired_m2,
       ST_NumGeometries(ST_MakeValid(g)) AS parts
FROM (VALUES (ST_GeomFromText('POLYGON((502000 2004000, 504000 2004000, 502000 2006000, 504000 2006000,
                                        502000 2004000))', 32643))) AS t (g);

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

The ring is invalid with the reason Self-intersection[1 1]. ST_MakeValid returns a MultiPolygon of both triangles, area 2; ST_Buffer(g, 0) returns one triangle, area 1. The manual mentions buffering by zero only as something “sometimes used” and points to ST_MakeValid as the better tool (ST_Buffer).

A repair keeps what the coordinates say, which is not always what was meant. Ward W02 typed with two corners swapped crosses itself at 503000 2005000; repaired, it is two triangles of 2,000,000 m² in all, not the 4 km² square of the ward. Check each repair against what you expect: the number of parts, the area, the ward’s neighbours.

Repairing by buffering with zero

UPDATE ward SET geom = ST_Buffer(geom, 0) turns the bow-tie into one triangle and reports nothing. Use ST_MakeValid(geom), then compare the areas before and after, so that a repair that lost or doubled a part shows up.

Checkpoint

Check: Finding and repairing invalid geometry

Question 1 of 2 What did ST_IsValidReason return for the bow-tie POLYGON((0 0, 2 2, 2 0, 0 2, 0 0))?
Choose one answer.

Ring Self-intersection[0 0] Valid Geometry Self-intersection[1 1]

Show the answer: Question 1 of 2

Answer: Self-intersection[1 1]

A ring that crosses itself is invalid, and ST_IsValidReason says where: at 1 1 for the bow-tie, at 503000 2005000 for the ward with two corners swapped.
Question 2 of 2 The bow-tie has two triangles of area 1 each. What did ST_Buffer(g, 0) return?
Choose one answer.

Both triangles, area 2 An empty polygon, area 0 One triangle, area 1, and no message

Show the answer: Question 2 of 2

Answer: One triangle, area 1, and no message

ST_MakeValid keeps every vertex and returned a MultiPolygon of area 2; the zero buffer returned one triangle. Compare areas before and after any repair.

Buffers and unions work in the units of the SRID

A buffer: The zone within a given distance of a feature, such as the area within 1 km of a clinic. Counting what lies inside a buffer answers many proximity questions. More in “What GIS is and the questions it answers” distance is in the units of the geometry’s spatial reference system, metres for the town’s UTM zone, and the circle is approximated with quad_segs “line segments used to approximate a quarter circle (default is 8)” (ST_Buffer). ST_Union merges many shapes into one, so an area two buffers share is counted once. The session buffers the 8 clinics by 1 km:

Buffers of 1 km around the clinics, their union, and what each covers

buffers_session.py

"""Buffers around the town's clinics, their union, and the households each method counts as covered."""
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("buffers.sql", db="town")

Output

$ 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 buffers.sql
-- One 1 km buffer: a polygon of 32 sides inside the circle, 8 segments for each quarter by default.
SELECT ST_NPoints(ST_Buffer(geom, 1000)) - 1 AS sides, round(ST_Area(ST_Buffer(geom, 1000))) AS buffer_m2,
       round(pi() * 1000 ^ 2) AS circle_m2
FROM clinic WHERE clinic = 'C2';
 sides | buffer_m2 | circle_m2 
-------+-----------+-----------
    32 |   3121445 |   3141593
(1 row)

-- The 8 clinics' buffers: their areas added up, and their union, which counts each overlap once.
SELECT round(sum(ST_Area(ST_Buffer(geom, 1000)))::numeric / 1e6, 3) AS sum_km2,
       round(ST_Area(ST_Union(ST_Buffer(geom, 1000)))::numeric / 1e6, 3) AS union_km2
FROM clinic;
 sum_km2 | union_km2 
---------+-----------
  24.972 |    24.587
(1 row)

-- Households covered: inside the union of buffers with 32 and 64 sides, and within 1,000 m by distance.
WITH u32 AS (SELECT ST_Union(ST_Buffer(geom, 1000)) AS g FROM clinic),
     u64 AS (SELECT ST_Union(ST_Buffer(geom, 1000, 'quad_segs=16')) AS g FROM clinic)
SELECT (SELECT count(*) FROM household h, u32 WHERE ST_Intersects(h.geom, u32.g)) AS in_32_sided,
       (SELECT count(*) FROM household h, u64 WHERE ST_Intersects(h.geom, u64.g)) AS in_64_sided,
       (SELECT count(*) FROM household h
        WHERE EXISTS (SELECT 1 FROM clinic c WHERE ST_DWithin(c.geom, h.geom, 1000))) AS within_1000_m;
 in_32_sided | in_64_sided | within_1000_m 
-------------+-------------+---------------
        2966 |        2971 |          2973
(1 row)

Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 buffers_session.py

buffers.sql

-- One 1 km buffer: a polygon of 32 sides inside the circle, 8 segments for each quarter by default.
SELECT ST_NPoints(ST_Buffer(geom, 1000)) - 1 AS sides, round(ST_Area(ST_Buffer(geom, 1000))) AS buffer_m2,
       round(pi() * 1000 ^ 2) AS circle_m2
FROM clinic WHERE clinic = 'C2';

-- The 8 clinics' buffers: their areas added up, and their union, which counts each overlap once.
SELECT round(sum(ST_Area(ST_Buffer(geom, 1000)))::numeric / 1e6, 3) AS sum_km2,
       round(ST_Area(ST_Union(ST_Buffer(geom, 1000)))::numeric / 1e6, 3) AS union_km2
FROM clinic;

-- Households covered: inside the union of buffers with 32 and 64 sides, and within 1,000 m by distance.
WITH u32 AS (SELECT ST_Union(ST_Buffer(geom, 1000)) AS g FROM clinic),
     u64 AS (SELECT ST_Union(ST_Buffer(geom, 1000, 'quad_segs=16')) AS g FROM clinic)
SELECT (SELECT count(*) FROM household h, u32 WHERE ST_Intersects(h.geom, u32.g)) AS in_32_sided,
       (SELECT count(*) FROM household h, u64 WHERE ST_Intersects(h.geom, u64.g)) AS in_64_sided,
       (SELECT count(*) FROM household h
        WHERE EXISTS (SELECT 1 FROM clinic c WHERE ST_DWithin(c.geom, h.geom, 1000))) AS within_1000_m;

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)")

Measure	Value
One buffer, 32 sides, against a true circle	3,121,445 m² against 3,141,593 m², 0.64 % short
The 8 buffers’ areas added up	24.972 km²
Their union	24.587 km², 0.385 km² of overlap counted once
Households inside the union of 32-sided buffers	2,966
Inside the union of 64-sided buffers	2,971
Within 1,000 m of a clinic, by `ST_DWithin`	2,973

A buffer is a polygon whose corners lie on the circle, so it misses the slivers between its sides and the arc: 7 households in reach of a clinic fell outside the 32-sided union.

Key idea

Count with ST_DWithin, which is exact; keep buffers for drawing and for areas, with more segments when the area matters.

The rest of the lesson goes further:

ST_Simplify drops vertices of each geometry on its own, and ST_SimplifyPreserveTopology does so while keeping each one valid.
ST_CoverageSimplify simplifies a set of polygons that share borders as one coverage (of polygons): A set of polygons that meet edge to edge, with no gaps and no overlaps, such as the districts of a state. Each shared border is stored in both neighbours, so tools that change borders must treat the set as a whole.. Neighbours then keep one common border.
ST_ClusterDBSCAN groups points by density, DBSCAN: Density-based clustering: points that lie close together, each with enough neighbours within a chosen distance, form clusters, and isolated points are left out as noise. The distance and the minimum count decide the result., with a distance and a minimum count. Points without enough neighbours stay out of every cluster.
ST_ClusterKMeans splits points into exactly k groups. You choose k, such as one group per clinic.
ST_SquareGrid and ST_HexagonGrid make cells to count in, and ST_GeoHash returns a geohash: A short text code for a cell of a world-wide grid of longitude and latitude. Each extra character names a smaller cell inside the previous one, so codes that share a beginning name nearby places., the name of a cell of a global grid. Cells of one size make the counts of different places comparable.

Checkpoint

Check: Buffers, unions and what they cover

Question 1 of 2 With the default quad_segs of 8, how many sides does a 1 km buffer of a point have?
Type a number.
Your answer sides

Show the answer: Question 1 of 2

Answer: 32 sides

8 segments for each quarter circle, 4 quarters: 32 sides. Its area was 3,121,445 m², 0.64 % short of the circle's 3,141,593 m².
Question 2 of 2 The union of 32-sided buffers held 2,966 households and ST_DWithin found 2,973 within 1,000 m. Which number belongs in a coverage report?
Choose one answer.

2,966, because the buffer is what the map shows Their average, 2,969.5 2,973, from ST_DWithin, which measures the real distance

Show the answer: Question 2 of 2

Answer: 2,973, from ST_DWithin, which measures the real distance

Count coverage with ST_DWithin; draw it with buffers, which may use more segments when the drawn area matters.

Simplifying shared borders without gaps or overlaps, measured Pro: opens with a Pro, Premium or Ultimate pass

Districts that share borders must be simplified together, or each side of a border moves on its own. This part builds a coverage of 30 districts, simplifies it three ways and measures the gaps and overlaps each way leaves.

In this part:

A coverage of 30 districts with staircase borders
Each district alone: the two copies of a border drift apart
The whole coverage at once: ST_CoverageSimplify

2 runnable examples (checks.sql, checks_session.py, coverage_session.py, districts.sql, simplify.sql) · about 6 minutes

Opening this part with your pass…

This part opens with any MySmartCoPilot pass.

Pro pass · 1 week₹79$2.99
30 days₹179$5.99
1 year₹1,699$45

One payment; a pass never renews. The price shown is the total you pay: no fees or charges are added at checkout.

Choose a pass I have a licence key

It also opens every other in-depth part in Learn and every Pro tool. Sales are final, except where the law gives you a right to withdraw, for example in the EU and the UK (refund policy).

Clusters, grids and geohashes on the town's households Pro: opens with a Pro, Premium or Ultimate pass

Points become information when they are grouped. This part clusters the 5,000 households by density and by k-means, counts them on square and hexagon grids, names cells with geohashes, and ends with where such processing belongs.

In this part:

DBSCAN: eps decides what a cluster is
k-means makes exactly k clusters, whatever the data
Grids and geohashes for counting and naming
Database, Python or QGIS?

2 runnable examples (clusters.sql, clusters_session.py, grids.sql, grids_session.py) · about 6 minutes

Opening this part with your pass…Pro part. It opens with the same pass as the part above.

Practise: find where a ring crosses itself

ST_IsValidReason told you where the bow-tie crossed itself. The exercise finds that point in Python, the way a validity check walks the edges of a ring:

Exercise · Medium · Python

Find the point where a ring crosses itself

ST_IsValidReason reported Self-intersection[1 1] for the bow-tie of this lesson. Write first_crossing(ring) in crossing.py: it returns the point (x, y) where two edges of the ring cross, or None when no two edges cross.

ring is a list of (x, y) points whose last point repeats the first, as in WKT.
The edges are the segments between consecutive points. Two edges that follow each other share a point, and so do the first and the last edge; those meetings are not crossings.
Test the pairs of edges in order, the first edge against every later one, then the second, and return the first crossing you find.

The sample tests use the bow-tie, the same bow-tie moved away from the origin, ward W02 typed with two corners swapped, and valid rings that must return None.

Starter code · crossing.py

def first_crossing(ring):
    """The first point (x, y) where two edges of a closed ring cross, or None."""
    # Replace this line: pair every edge with every later one that does not touch it at a shared point.
    return None

The sample tests · test_crossing.py

from crossing import first_crossing

BOWTIE = [(0, 0), (2, 2), (2, 0), (0, 2), (0, 0)]
W02_SWAPPED = [(502_000, 2_004_000), (504_000, 2_004_000), (502_000, 2_006_000), (504_000, 2_006_000), (502_000, 2_004_000)]


def test_the_bowtie_crosses_at_1_1():
    """finds the crossing that ST_IsValidReason reported for the bow-tie"""
    assert first_crossing(BOWTIE) == (1, 1)


def test_the_ward_with_two_corners_swapped():
    """finds the crossing of ward W02's two diagonals, at 503000 2005000"""
    assert first_crossing(W02_SWAPPED) == (503_000, 2_005_000)


def test_a_square_has_no_crossing():
    """returns None for a valid square, whose neighbouring edges only share their corners"""
    assert first_crossing([(0, 0), (2, 0), (2, 2), (0, 2), (0, 0)]) is None


def test_a_concave_ring_has_no_crossing():
    """returns None for a valid L-shaped ring with a reflex corner"""
    assert first_crossing([(0, 0), (3, 0), (3, 1), (1, 1), (1, 3), (0, 3), (0, 0)]) is None


def test_a_crossing_between_vertices():
    """finds a crossing that lies between the ring's vertices, at 2 1 for a stretched bow-tie"""
    assert first_crossing([(0, 0), (4, 2), (4, 0), (0, 2), (0, 0)]) == (2, 1)


def test_a_bowtie_away_from_the_origin():
    """finds the crossing of the same bow-tie moved one unit right and up, at 2 2"""
    assert first_crossing([(1, 1), (3, 3), (3, 1), (1, 3), (1, 1)]) == (2, 2)

A hint

For two segments p→q and r→s, write the points as p + t·(q − p) and r + u·(s − r). They cross when the two expressions are equal for some t and u between 0 and 1; solving the two equations gives t and u, with the denominator (q − p) × (s − r). When that denominator is 0 the segments are parallel and do not cross.

The sample tests run on this device, in your browser (Pyodide): nothing is sent to mysmartcopilot.com. The first run downloads Python (about 13.5 MB), which is kept for the next runs. A check in your browser is feedback for you, not proof that the code is right for every input.

Interview questions

Warm-up (fresher): Why is ST_Buffer(geom, 0) a risky way to repair an invalid polygon?

A buffer of zero rebuilds the polygon from what its edges enclose, and for a ring that crosses itself it can keep one part and drop another without any message. In the lesson the bow-tie, with two triangles, came back as one triangle with half the area. ST_MakeValid keeps every vertex and returned both triangles. Repair with it, and compare the area and the number of parts before and after.

Processing questions in interviews go further than repairs and buffers. These three start from this lesson, and the bank gives each one’s focus and opening move:

Mid to senior

A state's 30 district boundaries are simplified nightly with ST_SimplifyPreserveTopology for a dashboard. Every district is valid, yet analysts find slivers between neighbours and households counted in two districts. Explain why, and how you would produce simplified districts that still add up.

Model answer, Pro. Model answers open with any pass; prices are in the notice on this page. First step and follow-ups
Mid to senior

A city wants "hot spots" of dengue cases from 40,000 geocoded cases. DBSCAN with eps 50 m gives 2,000 tiny clusters; with eps 500 m, one cluster covers half the city. How do you choose eps and minpoints, and what do you report so that the result can be trusted?

Model answer, Pro. First step and follow-ups
Senior to staff

A team generalises the national boundary layer inside PostGIS every night for web maps, and analysts now complain that areas and counts changed. Which processing belongs in the database, and how do you keep analytical and display geometry apart without two teams editing two copies?

Model answer, Pro. First step and follow-ups

Flashcards · 6 cards

Recall it, then check

Say the answer to yourself before you look: recalling it is what makes it stick.

How do you find and repair invalid geometries in PostGIS?

ST_IsValidReason names the fault and its place, such as Self-intersection[1 1], and ST_MakeValid builds a valid geometry without losing any input vertex: a bow-tie came back as two triangles of area 2, where ST_Buffer(g, 0) kept one, area 1.
Why must a repaired geometry still be checked?

A repair keeps what the coordinates say, not what was meant: ward W02 typed with two corners swapped became two triangles of 2,000,000 m² in all, not its 4 km² square. Compare the parts, the area and the neighbours before and after.
In which units, and with how many segments, does ST_Buffer build a buffer?

In the units of the geometry's SRID, metres in the town's UTM zone, with quad_segs of 8 segments a quarter circle by default: a 1 km buffer of 32 sides covers 3,121,445 m² against a true circle's 3,141,593 m², 0.64 % short.
Why count coverage with ST_DWithin rather than with a union of buffers?

A buffer polygon misses the slivers between its sides and the arc: 2,966 households lay inside the union of 32-sided buffers and 2,971 inside the 64-sided one, while ST_DWithin found the exact 2,973 within 1,000 m. ST_Union counts shared areas once: 24.587 km² against 24.972 km² added up.
What do ST_SimplifyPreserveTopology, ST_CoverageSimplify, ST_ClusterDBSCAN and ST_ClusterKMeans do?

The first drops vertices of each geometry while keeping it valid; the second simplifies polygons that share borders as one coverage, so neighbours keep one border; DBSCAN groups points by density with a distance and a minimum count, leaving sparse points in no cluster; k-means splits points into exactly k groups.
Why count points in square, hexagonal or geohash cells?

Cells of one size, made by ST_SquareGrid or ST_HexagonGrid, or named by ST_GeoHash as cells of a global grid, make the counts of different places comparable.

Key takeaways

Find invalid rows with ST_IsValidReason, repair them with ST_MakeValid, and compare areas and parts before and after, because a repair keeps what was drawn, not what was meant.
Buffer in a projected SRID; a buffer is a 32-sided polygon by default, so count coverage with ST_DWithin and keep buffers for maps and areas.
ST_Union counts an area that buffers share once: 24.587 km² for the 8 buffers’ union, against 24.972 km² added up.

Terms in this lesson

coverage (of polygons): A set of polygons that meet edge to edge, with no gaps and no overlaps, such as the districts of a state. Each shared border is stored in both neighbours, so tools that change borders must treat the set as a whole.
DBSCAN: Density-based clustering: points that lie close together, each with enough neighbours within a chosen distance, form clusters, and isolated points are left out as noise. The distance and the minimum count decide the result.
geohash: A short text code for a cell of a world-wide grid of longitude and latitude. Each extra character names a smaller cell inside the previous one, so codes that share a beginning name nearby places.

References

PostGIS: ST_MakeValid (PostGIS Project Steering Committee)
PostGIS: ST_IsValidReason (PostGIS Project Steering Committee)
PostGIS: ST_Buffer (PostGIS Project Steering Committee)
PostGIS: ST_Union (PostGIS Project Steering Committee)
PostGIS: ST_SimplifyPreserveTopology (PostGIS Project Steering Committee)
PostGIS: ST_CoverageSimplify (PostGIS Project Steering Committee)
PostGIS: ST_CoverageInvalidEdges (PostGIS Project Steering Committee)
PostGIS: ST_ClusterDBSCAN (PostGIS Project Steering Committee)
PostGIS: ST_ClusterKMeans (PostGIS Project Steering Committee)
PostGIS: ST_SquareGrid (PostGIS Project Steering Committee)
PostGIS: ST_GeoHash (PostGIS Project Steering Committee)
GEOS: CoverageSimplifier (GEOS project (OSGeo))

Report a problem with this lesson

In this module: Spatial SQL with PostGIS

PostGIS setup: geometry, geography and SRIDs Free lesson, Pro answers: 15 minutes
Spatial relationships and joins in SQL Free core, Pro depth: 6 of 15 minutes free
PostGIS indexes, KNN and query plans Free core, Pro depth: 6 of 15 minutes free
Geometry processing in PostGIS Free core, Pro depth: 6 of 15 minutes free
Serving GeoJSON, vector tiles and OGC APIs Free core, Pro depth: 6 of 15 minutes free

Your country

Geometry processing in PostGIS

What you will learn

Before you start

What you will learn

Before you start

Repair first: name the fault, then keep what was drawn

Buffers and unions work in the units of the SRID

Simplifying shared borders without gaps or overlaps, measured Pro: opens with a Pro, Premium or Ultimate pass

Clusters, grids and geohashes on the town's households Pro: opens with a Pro, Premium or Ultimate pass

Practise: find where a ring crosses itself

Find the point where a ring crosses itself

Results of the sample tests

Interview questions

A state's 30 district boundaries are simplified nightly with ST_SimplifyPreserveTopology for a dashboard. Every district is valid, yet analysts find slivers between neighbours and households counted in two districts. Explain why, and how you would produce simplified districts that still add up.

A city wants "hot spots" of dengue cases from 40,000 geocoded cases. DBSCAN with eps 50 m gives 2,000 tiny clusters; with eps 500 m, one cluster covers half the city. How do you choose eps and minpoints, and what do you report so that the result can be trusted?

A team generalises the national boundary layer inside PostGIS every night for web maps, and analysts now complain that areas and counts changed. Which processing belongs in the database, and how do you keep analytical and display geometry apart without two teams editing two copies?

Key takeaways

References

Related tools