feat: add 'urban_centres' dataset for better city coverage and add population ranking

DEPLOY_TO_HF.md

@@ -91,11 +91,13 @@ Ensure normalized data exists locally:
 data/overture_normalized/divisions_area/*.parquet
 data/overture_normalized/localities/*.parquet
 data/natural_earth_normalized/ne_geography.parquet
+data/urban_centres_normalized/*.parquet
 ```
 
-`localities/` holds cities and towns. Only fuzzy search reads it; the
-natural-language pipeline reads `divisions_area/` alone, because the model
-was not trained on those subtypes.
+`localities/` holds cities and towns, and `urban_centres_normalized/` the
+GHSL urban centres. Only fuzzy search reads them; the natural-language
+pipeline reads `divisions_area/` alone, because the model was not trained on
+those subtypes.
 
 ### Upload dataset
 
@@ -113,6 +115,13 @@ hf upload developmentseed/gazet-geodata \
     /natural_earth_normalized \
     --repo-type dataset \
     --commit-message <message>
+
+# Upload GHSL urban centres
+hf upload developmentseed/gazet-geodata \
+    ./data/urban_centres_normalized \
+    /urban_centres_normalized \
+    --repo-type dataset \
+    --commit-message <message>
 ```
 
 The HF repo will only contain:
@@ -124,8 +133,10 @@ gazet-geodata/
 │   │   └── *.parquet
 │   └── localities/
 │       └── *.parquet
-└── natural_earth_normalized/
-    └── ne_geography.parquet
+├── natural_earth_normalized/
+│   └── ne_geography.parquet
+└── urban_centres_normalized/
+    └── *.parquet
 ```
 
 This keeps the repo size minimal - only what's needed for deployment.
@@ -218,6 +229,12 @@ hf upload developmentseed/gazet-geodata \
     /natural_earth_normalized \
     --repo-type dataset \
     --commit-message <message>
+
+hf upload developmentseed/gazet-geodata \
+    ./data/urban_centres_normalized \
+    /urban_centres_normalized \
+    --repo-type dataset \
+    --commit-message <message>
 ```
 
 The Space will use the updated data on next rebuild.

README.md

@@ -31,15 +31,18 @@ uv sync --extra dev --extra demo
 
 ### Data preparation
 
-1. Download Overture divisions data
+1. Download Overture divisions data: the areas, and the division records that carry each place's population
 2. Download the 10m physical layer from [Natural Earth](https://www.naturalearthdata.com/downloads/10m-physical-vectors/)
 3. Unzip the data
 4. Convert natural earth data to parquet
+5. Download the [GHSL Urban Centre Database](https://human-settlement.emergency.copernicus.eu/ghs_ucdb_2024.php) (GHS-UCDB R2024A): named city polygons for places Overture has only as a point, such as London
+6. Normalize everything: `python -m dataset.scripts.normalize_geodata`
 
 Example for downloading overture
 
 ```bash
 aws s3 sync s3://overturemaps-us-west-2/release/2026-02-18.0/theme=divisions/type=division_area/ data/overture/divisions_area
+aws s3 sync s3://overturemaps-us-west-2/release/2026-02-18.0/theme=divisions/type=division/ data/overture/division
 ```
 
 Example for running conversion script for natural earth
@@ -49,6 +52,13 @@ unzip ~/Downloads/10m_physical.zip -d data/natural_earth
 python -m ingest.convert_natural_earth data/natural_earth
 ```
 
+Example for downloading the GHSL urban centres
+
+```bash
+curl -LO https://jeodpp.jrc.ec.europa.eu/ftp/jrc-opendata/GHSL/GHS_UCDB_GLOBE_R2024A/GHS_UCDB_GLOBE_R2024A/V1-1/GHS_UCDB_GLOBE_R2024A_V1_1.zip
+unzip GHS_UCDB_GLOBE_R2024A_V1_1.zip GHS_UCDB_GLOBE_R2024A.gpkg -d data/ghsl
+```
+
 ### Based on ollama
 
 For now, gazet relies on [ollama](https://ollama.com/). For remote (cloud) models, ensure you are loged into Ollama.
@@ -90,8 +100,9 @@ uv run streamlit run gazet_demo.py   # demo UI
 ## Design notes
 
 - `api.py` exposes GET `/search?q=<query>`; returns GeoJSON FeatureCollection and logs intermediate output. This path runs two LLM calls (place extraction, then SQL generation for candidate selection/disambiguation and GIS ops).
-- `api.py` also exposes GET `/search/fuzzy?q=<place name>&limit=5&simplify=true&sources=divisions_area,natural_earth` — pure Jaro-Winkler fuzzy match with no LLM involved. `q` is a place-name string, not a natural-language query (there's no place-extraction step). Returns the combined top-`limit` matches across sources as a GeoJSON FeatureCollection with raw (or simplified, if `simplify=true`) geometry attached. Pass `ids_only=true` to skip full geometry and get back `{"ids": [{"source", "id", "name", "matched_name", "country", "subtype", "admin_level", "bbox", "similarity", "is_substring_match"}, ...]}` — `country`/`subtype`/`admin_level` disambiguate same-named places (e.g. multiple real "Loja"s across Ecuador and Spain, or Ecuador's "Loja" region vs. its nested "Loja" county), `bbox` is `[minx, miny, maxx, maxy]`, a much smaller payload than full geometry, for picking a candidate before fetching its full geometry, and `similarity` is the score the list was ranked on — Jaro-Winkler always returns a best row, so without it a caller cannot tell a hit from the closest thing to a miss.
-- Fuzzy search (`mode=fuzzy`) also searches cities and towns (Overture `locality` and `localadmin`), stored apart in `overture_normalized/localities/`. The natural-language mode never sees them, because its model was not trained on them.
+- `api.py` also exposes GET `/search/fuzzy?q=<place name>&limit=5&simplify=true&sources=divisions_area,natural_earth,urban_centres` — pure Jaro-Winkler fuzzy match with no LLM involved. `q` is a place-name string, not a natural-language query (there's no place-extraction step). Returns the combined top-`limit` matches across sources as a GeoJSON FeatureCollection with raw (or simplified, if `simplify=true`) geometry attached. Pass `ids_only=true` to skip full geometry and get back `{"ids": [{"source", "id", "name", "matched_name", "country", "subtype", "admin_level", "bbox", "similarity", "is_substring_match"}, ...]}` — `country`/`subtype`/`admin_level` disambiguate same-named places (e.g. multiple real "Loja"s across Ecuador and Spain, or Ecuador's "Loja" region vs. its nested "Loja" county), `bbox` is `[minx, miny, maxx, maxy]`, a much smaller payload than full geometry, for picking a candidate before fetching its full geometry, and `similarity` is the score the list was ranked on — Jaro-Winkler always returns a best row, so without it a caller cannot tell a hit from the closest thing to a miss.
+- Fuzzy search (`mode=fuzzy`) also searches cities and towns (Overture `locality` and `localadmin`), stored apart in `overture_normalized/localities/`, and GHSL urban centres (`urban_centres_normalized/`, ids prefixed `ghsl_`): named city polygons for places Overture has only as a point, such as London. The natural-language mode never sees either, because its model was not trained on them.
+- Among equal scores, fuzzy search puts the more populous place first, so "London" returns the UK's before a same-named town. Population comes from Overture's `division` records and from GHSL; a place without one sorts last.
 - Fuzzy matching (`search.py`) ranks exact-substring matches (query literally contained in the name, e.g. "Loja" in "Loja Province") ahead of pure edit-distance near-misses (e.g. "Lodja"), even when the near-miss scores a slightly higher raw Jaro-Winkler similarity. Each record is matched under two names and keeps the better score: its own (`names.primary`) and its English one (`names.common.en` for `divisions_area`, `names.en` for `natural_earth`). Neither works alone — `names.common.en` is null for 77% of `divisions_area` rows, which matching on English alone would drop, while matching on `names.primary` alone leaves a place unreachable by the English name it is usually asked for ("Copenhagen" for `Københavns Kommune`). `matched_name` reports which of the two the query hit; `name` is the English one where the record has one, the same expression `/geometry/{id}` uses, so a candidate and a later lookup of it agree on what the place is called.
 - `api.py` also exposes GET `/geometry/{id}?source=divisions_area&simplify=true` — fetch a single feature's geometry directly by ID, no fuzzy matching or LLM involved. `source` is inferred from the ID if omitted (Natural Earth IDs are prefixed `ne_`). Returns a single GeoJSON Feature, 404 if not found.
 - The spatial extension is loaded once at API startup (FastAPI `lifespan`) onto a shared DuckDB connection stored on `app.state`; each request gets a cheap `cursor()` off it instead of paying DuckDB's ~90ms `LOAD spatial` cost per call.

dataset/scripts/normalize_geodata.py

@@ -13,6 +13,7 @@ Output layout under data/ by default:
     overture_normalized/divisions_area/part-000.parquet
     overture_normalized/localities/part-000.parquet
     natural_earth_normalized/ne_geography.parquet
+    urban_centres_normalized/part-000.parquet
 
 Cities and towns (Overture ``locality`` and ``localadmin``) go to their own
 file. Fuzzy search reads it; the training pipeline and the natural-language
@@ -21,7 +22,14 @@ those subtypes.
 
 Both Overture files also store ``search_names``: every name a record
 answers to, folded for comparison, so fuzzy search does not rebuild and
-fold them on every query.
+fold them on every query. They also store ``population``, joined from
+Overture's ``division`` records, so fuzzy search can rank a city above a
+same-named village.
+
+Urban centres come from the GHSL Urban Centre Database (GHS-UCDB R2024A):
+named city polygons for places Overture has only as a point, such as
+London. They are a source of their own, never merged into an Overture
+record, so every result traces back to one dataset.
 """
 
 from pathlib import Path
@@ -29,12 +37,24 @@ from pathlib import Path
 import duckdb
 
 from gazet.config import _DATA_DIR
-from gazet.search import DIVISIONS_AREA_SEARCH_NAMES, SEARCH_NAMES_COLUMN
+from gazet.search import (
+    DIVISIONS_AREA_SEARCH_NAMES,
+    SEARCH_NAMES_COLUMN,
+    search_names_sql,
+)
 
 #: The subtypes the natural-language model was trained on.
 TRAINED_SUBTYPES = ("country", "region", "county")
 #: Cities and towns, for fuzzy search only.
 LOCALITY_SUBTYPES = ("localadmin", "locality")
+#: The GHS-UCDB layer that carries each urban centre's names, country and
+#: population. The other layers repeat the same polygons with other attributes.
+GHSL_UCDB_LAYER = "GHS_UCDB_THEME_GENERAL_CHARACTERISTICS_GLOBE_R2024A"
+
+
+def _ghsl(column: str) -> str:
+    """A GHS-UCDB column, whose name and text values start with a byte-order mark."""
+    return f'ltrim("\ufeff{column}", chr(65279))'
 
 
 def normalize_geodata(output_root: Path | None = None) -> dict[str, str]:
@@ -51,7 +71,13 @@ def normalize_geodata(output_root: Path | None = None) -> dict[str, str]:
     overture_dir = root / "overture_normalized" / "divisions_area"
     localities_dir = root / "overture_normalized" / "localities"
     natural_earth_dir = root / "natural_earth_normalized"
-    for directory in (overture_dir, localities_dir, natural_earth_dir):
+    urban_centres_dir = root / "urban_centres_normalized"
+    for directory in (
+        overture_dir,
+        localities_dir,
+        natural_earth_dir,
+        urban_centres_dir,
+    ):
         directory.mkdir(parents=True, exist_ok=True)
 
     natural_earth_path = natural_earth_dir / "ne_geography.parquet"
@@ -69,9 +95,14 @@ def normalize_geodata(output_root: Path | None = None) -> dict[str, str]:
             f"""
             COPY (
                 SELECT
-                    * REPLACE (ST_GeomFromWKB(ST_AsWKB(geometry)) AS geometry),
+                    area.* REPLACE (ST_GeomFromWKB(ST_AsWKB(area.geometry)) AS geometry),
+                    division.population,
                     {DIVISIONS_AREA_SEARCH_NAMES} AS {SEARCH_NAMES_COLUMN}
-                FROM read_parquet('{root / "overture/divisions_area/*.parquet"}')
+                FROM read_parquet('{root / "overture/divisions_area/*.parquet"}') AS area
+                LEFT JOIN (
+                    SELECT id, population
+                    FROM read_parquet('{root / "overture/division/*.parquet"}')
+                ) AS division ON division.id = area.division_id
                 WHERE geometry IS NOT NULL
                   AND subtype IN {subtypes}
                   AND is_land = true
@@ -90,12 +121,57 @@ def normalize_geodata(output_root: Path | None = None) -> dict[str, str]:
         ) TO '{natural_earth_path}' (FORMAT PARQUET)
         """
     )
+
+    # Shaped like the Overture files, so search and lookup read every source
+    # the same way. Names come as one "; "-separated list of alternates.
+    alternates = f"string_split({_ghsl('GC_UCN_LIS_2025')}, '; ')"
+    con.execute(
+        f"""
+        COPY (
+            SELECT
+                * EXCLUDE (alternate_names),
+                {search_names_sql("names.primary", "CAST(NULL AS VARCHAR)", "alternate_names")}
+                    AS {SEARCH_NAMES_COLUMN}
+            FROM (
+                SELECT
+                    'ghsl_' || CAST("\ufeffID_UC_G0" AS VARCHAR) AS id,
+                    ST_GeomFromWKB(ST_AsWKB(geometry)) AS geometry,
+                    {{
+                        'xmin': ST_XMin(geometry), 'xmax': ST_XMax(geometry),
+                        'ymin': ST_YMin(geometry), 'ymax': ST_YMax(geometry)
+                    }} AS bbox,
+                    {_ghsl("GC_CNT_GAD_2025")} AS country,
+                    'urban_centre' AS subtype,
+                    CAST(NULL AS VARCHAR) AS class,
+                    {{'primary': {_ghsl("GC_UCN_MAI_2025")}}} AS names,
+                    CAST(NULL AS VARCHAR) AS region,
+                    CAST(NULL AS INTEGER) AS admin_level,
+                    true AS is_land,
+                    CAST(NULL AS BOOLEAN) AS is_territorial,
+                    CAST(round("\ufeffGC_POP_TOT_2025") AS BIGINT) AS population,
+                    {alternates} AS alternate_names
+                FROM (
+                    SELECT
+                        *,
+                        -- Mollweide to lon/lat.
+                        ST_Transform(geom, 'ESRI:54009', 'EPSG:4326', always_xy := true)
+                            AS geometry
+                    FROM st_read(
+                        '{root / "ghsl/GHS_UCDB_GLOBE_R2024A.gpkg"}',
+                        layer = '{GHSL_UCDB_LAYER}'
+                    )
+                )
+            )
+        ) TO '{urban_centres_dir / "part-000.parquet"}' (FORMAT PARQUET)
+        """
+    )
     con.close()
 
     return {
         "divisions_area": str(overture_dir / "*.parquet"),
         "localities": str(localities_dir / "*.parquet"),
         "natural_earth": str(natural_earth_path),
+        "urban_centres": str(urban_centres_dir / "*.parquet"),
     }
 
 

src/gazet/api.py

@@ -32,7 +32,8 @@ from .sql import run_geo_sql_dspy, run_geo_sql_gguf
 #: endpoint reports the backends, not the code.
 API_VERSION = version("gazet")
 
-_FUZZY_SOURCES = ("divisions_area", "natural_earth")
+#: urban_centres is fuzzy-only: the natural-language model was not trained on it.
+_FUZZY_SOURCES = ("divisions_area", "natural_earth", "urban_centres")
 
 SearchMode = Literal["nl", "fuzzy"]
 
@@ -61,8 +62,8 @@ _Q_SIMPLIFY = Query(
 _Q_SOURCES = Query(
     None,
     description=(
-        "Comma-separated subset of divisions_area,natural_earth for mode=fuzzy "
-        "(defaults to both); ignored for mode=nl."
+        "Comma-separated subset of divisions_area,natural_earth,urban_centres "
+        "for mode=fuzzy (defaults to all); ignored for mode=nl."
     ),
 )
 _Q_IDS_ONLY = Query(
@@ -293,13 +294,16 @@ def health(request: Request) -> dict[str, Any]:
 def sources(request: Request) -> dict[str, Any]:
     """List available data sources with row counts and name ranges."""
     con = request.app.state.duckdb_con
-    from .config import DIVISIONS_AREA_PATH, NATURAL_EARTH_PATH
+    from .config import DIVISIONS_AREA_PATH, NATURAL_EARTH_PATH, URBAN_CENTRES_PATH
 
     info = {}
     for name, path in [
         ("divisions_area", DIVISIONS_AREA_PATH),
         ("natural_earth", NATURAL_EARTH_PATH),
+        ("urban_centres", URBAN_CENTRES_PATH),
     ]:
+        if not path:
+            continue
         try:
             row = con.execute(
                 f"SELECT COUNT(*) as count, MIN(names.primary) as min_name, MAX(names.primary) as max_name FROM read_parquet('{path}')"
@@ -362,9 +366,10 @@ def _fuzzy_search(
     ``q`` is a place-name string (not a natural-language query) — no
     place-extraction step. Matches are ranked by Jaro-Winkler similarity
     across the requested ``sources`` (comma-separated subset of
-    divisions_area/natural_earth; defaults to both), combined and truncated
-    to the top ``limit``. Each record is matched under both its own name and
-    its English one, so "Copenhagen" and "København" find the same place.
+    divisions_area/natural_earth/urban_centres; defaults to all), combined
+    and truncated to the top ``limit``, the more populous place first among
+    equal scores. Each record is matched under both its own name and its
+    English one, so "Copenhagen" and "København" find the same place.
 
     Pass ``ids_only=true`` to skip fetching full geometry and get back
     lightweight candidates (id/name/country/subtype/admin_level/bbox, plus
@@ -403,7 +408,13 @@ def _fuzzy_search(
         candidates_df = (
             pd.concat(candidate_dfs, ignore_index=True)
             .drop_duplicates(subset=["source", "id"])
-            .sort_values(["is_substring_match", "similarity"], ascending=[False, False])
+            # Stable, so rows that tie on all three keep each source's order.
+            .sort_values(
+                ["is_substring_match", "similarity", "population"],
+                ascending=[False, False, False],
+                na_position="last",
+                kind="stable",
+            )
             .head(limit)
             .reset_index(drop=True)
         )
@@ -533,7 +544,7 @@ def search_fuzzy(
     simplify: bool = Query(True, description="Simplify geometry to GeoJSON."),
     sources: str | None = Query(
         None,
-        description="Comma-separated subset of divisions_area,natural_earth (defaults to both).",
+        description="Comma-separated subset of divisions_area,natural_earth,urban_centres (defaults to all).",
     ),
     ids_only: bool = Query(
         False,
@@ -558,16 +569,17 @@ def get_geometry(
     id: str,
     source: str | None = Query(
         None,
-        description="Restrict lookup to 'divisions_area' or 'natural_earth'; inferred from id if omitted.",
+        description="Restrict lookup to 'divisions_area', 'natural_earth' or 'urban_centres'; inferred from id if omitted.",
     ),
     simplify: bool = Query(True, description="Simplify geometry to GeoJSON."),
 ) -> dict[str, Any]:
     """Fetch a single feature's geometry directly by ID — no fuzzy matching, no LLM.
 
-    ``source`` restricts the lookup to ``divisions_area`` or ``natural_earth``;
-    if omitted, it's inferred from the ID (Natural Earth IDs are prefixed
-    ``ne_``). Returns a single GeoJSON Feature, or 404 if the ID doesn't exist
-    in the resolved source.
+    ``source`` restricts the lookup to ``divisions_area``, ``natural_earth``
+    or ``urban_centres``; if omitted, it's inferred from the ID (Natural
+    Earth IDs are prefixed ``ne_``, urban centre IDs ``ghsl_``). Returns a
+    single GeoJSON Feature, or 404 if the ID doesn't exist in the resolved
+    source.
     """
     if source is not None and source not in _FUZZY_SOURCES:
         raise HTTPException(status_code=400, detail=f"Unknown source: {source}")

src/gazet/config.py

@@ -32,16 +32,24 @@ DIVISIONS_AREA_PATH = str(
         _DATA_DIR / "overture/divisions_area/*.parquet",
     )
 )
+
+
+def _normalized_only(glob: pathlib.Path) -> str:
+    """A dataset that exists only as a normalized copy; empty where it was never built."""
+    if os.environ.get("GAZET_USE_NORMALIZED_DATA", "1") == "0":
+        return ""
+    return str(glob) if any(glob.parent.glob(glob.name)) else ""
+
+
 # Cities and towns (Overture locality and localadmin), kept out of
 # DIVISIONS_AREA_PATH because the natural-language model was not trained on
-# them. Only fuzzy search reads them; empty where they were never built.
-_LOCALITIES_GLOB = _DATA_DIR / "overture_normalized/localities/*.parquet"
-LOCALITIES_PATH = (
-    str(_LOCALITIES_GLOB)
-    if os.environ.get("GAZET_USE_NORMALIZED_DATA", "1") != "0"
-    and any(_LOCALITIES_GLOB.parent.glob(_LOCALITIES_GLOB.name))
-    else ""
+# them. Only fuzzy search reads them.
+LOCALITIES_PATH = _normalized_only(
+    _DATA_DIR / "overture_normalized/localities/*.parquet"
 )
+# GHSL urban centres: named city polygons for places Overture has only as a
+# point, such as London. Only fuzzy search reads them.
+URBAN_CENTRES_PATH = _normalized_only(_DATA_DIR / "urban_centres_normalized/*.parquet")
 NATURAL_EARTH_PATH = str(
     _prefer_normalized(
         _DATA_DIR / "natural_earth_normalized/ne_geography.parquet",

src/gazet/search.py

@@ -4,7 +4,12 @@ from functools import cache, partial
 import duckdb
 import pandas as pd
 
-from .config import DIVISIONS_AREA_PATH, LOCALITIES_PATH, NATURAL_EARTH_PATH
+from .config import (
+    DIVISIONS_AREA_PATH,
+    LOCALITIES_PATH,
+    NATURAL_EARTH_PATH,
+    URBAN_CENTRES_PATH,
+)
 from .schemas import Place
 
 logger = logging.getLogger(__name__)
@@ -73,24 +78,37 @@ DIVISIONS_AREA_SEARCH_NAMES = search_names_sql(
 NATURAL_EARTH_SEARCH_NAMES = search_names_sql(PRIMARY_NAME, NATURAL_EARTH_ENGLISH_NAME)
 
 
+#: The column normalize_geodata stores a place's population in, where its
+#: source has one: Overture joined from its ``division`` records, GHSL from
+#: its own. Ranks a city above a same-named village.
+POPULATION_COLUMN = "population"
+NO_POPULATION = "CAST(NULL AS BIGINT)"
+
+
 @cache
-def _stores_search_names(path: str) -> bool:
-    """Whether a parquet file (or glob) carries the stored search names."""
+def _stored_columns(path: str) -> frozenset[str]:
+    """The columns a parquet file (or glob) carries."""
     columns = duckdb.connect().execute("DESCRIBE SELECT * FROM read_parquet(?)", [path])
-    return SEARCH_NAMES_COLUMN in {row[0] for row in columns.fetchall()}
+    return frozenset(row[0] for row in columns.fetchall())
 
 
-def _search_names(path: str | list[str], computed: str) -> str:
-    """The stored search names where every file has them, else ``computed``.
+def _stored(path: str | list[str], column: str, computed: str) -> str:
+    """``column`` where every file stores it, else ``computed``.
 
-    Raw Overture downloads and the test fixtures have no stored column.
+    Raw Overture downloads, files normalized before a column was added, and
+    the test fixtures may lack one.
     """
     paths = [path] if isinstance(path, str) else path
-    if all(_stores_search_names(p) for p in paths):
-        return SEARCH_NAMES_COLUMN
+    if all(column in _stored_columns(p) for p in paths):
+        return column
     return computed
 
 
+def _search_names(path: str | list[str], computed: str) -> str:
+    """The stored search names where every file has them, else ``computed``."""
+    return _stored(path, SEARCH_NAMES_COLUMN, computed)
+
+
 def _divisions_area_paths(include_localities: bool) -> str | list[str]:
     """The divisions_area file, and the localities beside it where asked for.
 
@@ -115,6 +133,7 @@ def simple_fuzzy_search(
     include_geometry: bool = False,
     include_bbox: bool = False,
     bbox_expr: str = GEOMETRY_BBOX,
+    population_expr: str = NO_POPULATION,
 ) -> pd.DataFrame:
     """Jaro-Winkler similarity search using only the place name.
 
@@ -137,6 +156,10 @@ def simple_fuzzy_search(
     Jaro-Winkler similarity — substring containment is a stronger signal of
     a true match than character-level similarity alone.
 
+    Among equal scores, the larger ``population_expr`` comes first, so
+    "London" puts the city of ten million above a village of the same name.
+    A place with no population sorts after every place that has one.
+
     ``include_geometry``/``include_bbox`` require the spatial extension to
     already be loaded on ``con`` (``INSTALL spatial; LOAD spatial;``).
     ``include_bbox`` returns ``[minx, miny, maxx, maxy]`` from ``bbox_expr``
@@ -163,6 +186,7 @@ def simple_fuzzy_search(
                 admin_level,
                 is_land,
                 is_territorial{extra_clause}{bbox_clause}{carried_geometry},
+                CAST({population_expr} AS BIGINT) AS population,
                 {_readable(name_expr)} AS primary_name,
                 {_readable(english_expr)} AS english_name,
                 {search_names_expr} AS search_names
@@ -192,7 +216,8 @@ def simple_fuzzy_search(
             -- id last, so rows that tie on everything else (such as
             -- same-named towns, which have no admin_level) come back
             -- in the same order on every run.
-            ORDER BY best.substring DESC, best.similarity DESC, admin_level ASC, id
+            ORDER BY best.substring DESC, best.similarity DESC,
+                population DESC NULLS LAST, admin_level ASC, id
             LIMIT $3
         )
         SELECT
@@ -205,11 +230,13 @@ def simple_fuzzy_search(
             admin_level,
             is_land,
             is_territorial{extra_clause}{", bbox" if include_bbox else ""}{geometry_clause},
+            population,
             best.name AS matched_name,
             best.similarity AS similarity,
             best.substring AS is_substring_match
         FROM ranked
-        ORDER BY is_substring_match DESC, similarity DESC, admin_level ASC, id
+        ORDER BY is_substring_match DESC, similarity DESC,
+            population DESC NULLS LAST, admin_level ASC, id
         """,
         [place.place, path, limit],
     )
@@ -254,6 +281,35 @@ def search_divisions_area(
         include_geometry=include_geometry,
         include_bbox=include_bbox,
         bbox_expr=DIVISIONS_AREA_BBOX,
+        population_expr=_stored(path, POPULATION_COLUMN, NO_POPULATION),
+    )
+
+
+def search_urban_centres(
+    con: duckdb.DuckDBPyConnection,
+    place: Place,
+    limit: int = 5,
+    include_geometry: bool = False,
+    include_bbox: bool = False,
+) -> pd.DataFrame:
+    """Fuzzy-match a place against GHSL urban centres (named city polygons).
+
+    Matched on each centre's main name and its alternates. Returns an empty
+    frame where the urban centres were never built.
+    """
+    if not URBAN_CENTRES_PATH:
+        return pd.DataFrame()
+    return simple_fuzzy_search(
+        con,
+        URBAN_CENTRES_PATH,
+        "urban_centres",
+        place,
+        SEARCH_NAMES_COLUMN,
+        limit=limit,
+        include_geometry=include_geometry,
+        include_bbox=include_bbox,
+        bbox_expr=DIVISIONS_AREA_BBOX,
+        population_expr=POPULATION_COLUMN,
     )
 
 
@@ -356,16 +412,32 @@ def get_natural_earth_by_id(
     )
 
 
+def get_urban_centre_by_id(
+    con: duckdb.DuckDBPyConnection, id: str, include_geometry: bool = True
+) -> pd.DataFrame:
+    """Look up a single GHSL urban centre by exact ID; empty where none were built."""
+    if not URBAN_CENTRES_PATH:
+        return pd.DataFrame()
+    return fetch_by_id(
+        con, URBAN_CENTRES_PATH, "urban_centres", id, include_geometry=include_geometry
+    )
+
+
 _SOURCE_SEARCH_FNS = {
     "divisions_area": search_divisions_area,
     "natural_earth": search_natural_earth,
+    "urban_centres": search_urban_centres,
 }
 
 _SOURCE_FETCH_FNS = {
     "divisions_area": get_division_by_id,
     "natural_earth": get_natural_earth_by_id,
+    "urban_centres": get_urban_centre_by_id,
 }
 
+#: ID prefixes that name their source; any other ID is a divisions_area one.
+_ID_PREFIXES = {"ne_": "natural_earth", "ghsl_": "urban_centres"}
+
 
 def get_by_id(
     con: duckdb.DuckDBPyConnection,
@@ -375,11 +447,13 @@ def get_by_id(
 ) -> pd.DataFrame:
     """Look up a single record by ID, inferring the source if not given.
 
-    Natural Earth IDs are always prefixed ``ne_`` (see ``config.SCHEMA_INFO``);
-    anything else is assumed to be a divisions_area ID.
+    Natural Earth IDs are always prefixed ``ne_`` (see ``config.SCHEMA_INFO``)
+    and urban centre IDs ``ghsl_``; anything else is assumed to be a
+    divisions_area ID.
     """
-    resolved_source = source or (
-        "natural_earth" if id.startswith("ne_") else "divisions_area"
+    resolved_source = source or next(
+        (src for prefix, src in _ID_PREFIXES.items() if id.startswith(prefix)),
+        "divisions_area",
     )
     return _SOURCE_FETCH_FNS[resolved_source](
         con, id, include_geometry=include_geometry

tests/conftest.py

@@ -87,8 +87,27 @@ _LOCALITY_ROWS = [
     # reverse id order, so a search only returns them sorted if it sorts.
     ("lisbon-b", "Lisbon", None, "US", "locality", None, -70.1, 44.0),
     ("lisbon-a", "Lisbon", None, "US", "locality", None, -83.2, 42.9),
+    # Same name, but populations tell them apart: the larger comes first,
+    # whatever its id.
+    ("london-ky", "London", None, "US", "locality", None, -84.1, 37.1),
+    ("london-oh", "London", None, "US", "locality", None, -83.4, 39.9),
 ]
 
+#: GHSL urban centres, shaped like the normalized urban_centres file. The
+#: UK's London is the place Overture has only as a point.
+_URBAN_CENTRE_ROWS = [
+    ("ghsl_2255", "London", None, "Canada", "urban_centre", None, -81.3, 42.9),
+    ("ghsl_5816", "London", None, "United Kingdom", "urban_centre", None, -0.5, 51.3),
+]
+
+#: ``population`` where the fixture has one; every other row has none.
+_POPULATION = {
+    "london-ky": 7_000,
+    "london-oh": 10_000,
+    "ghsl_2255": 365_269,
+    "ghsl_5816": 10_408_332,
+}
+
 
 def _literal(value):
     """A string as a SQL literal, or NULL where there is no value."""
@@ -121,7 +140,8 @@ def _write_divisions(path, rows):
     values = ",\n".join(
         f"({_literal(id)}, {_literal(primary)}, {_name_map(english)}, {_rules(id)}, "
         f"{_literal(country)}, {_literal(subtype)}, "
-        f"{'NULL' if admin_level is None else admin_level}, {lon}, {lat})"
+        f"{'NULL' if admin_level is None else admin_level}, {lon}, {lat}, "
+        f"{_POPULATION.get(id, 'NULL')})"
         for id, primary, english, country, subtype, admin_level, lon, lat in rows
     )
     c.execute(
@@ -149,11 +169,13 @@ def _write_divisions(path, rows):
                 CAST(admin_level AS INTEGER) AS admin_level,
                 true AS is_land,
                 false AS is_territorial,
-                id AS division_id
+                id AS division_id,
+                CAST(population AS BIGINT) AS population
             FROM (VALUES
                 {values}
             ) AS t(
-                id, primary_name, common, rules, country, subtype, admin_level, lon, lat
+                id, primary_name, common, rules, country, subtype, admin_level,
+                lon, lat, population
             )
             )
         ) TO '{path}' (FORMAT PARQUET)
@@ -177,6 +199,20 @@ def localities_parquet(tmp_path_factory):
     return _write_divisions(path, _LOCALITY_ROWS)
 
 
+@pytest.fixture(scope="session")
+def urban_centres_parquet(tmp_path_factory):
+    """Write the urban centre rows to a parquet with stored search names."""
+    path = tmp_path_factory.mktemp("urban_centres") / "urban_centres.parquet"
+    return _write_divisions(path, _URBAN_CENTRE_ROWS)
+
+
+@pytest.fixture()
+def urban_centres_source(monkeypatch, urban_centres_parquet):
+    """Point urban_centres at the fixture parquet for the duration of a test."""
+    monkeypatch.setattr("gazet.search.URBAN_CENTRES_PATH", urban_centres_parquet)
+    return urban_centres_parquet
+
+
 @pytest.fixture()
 def exonym_source(monkeypatch, exonym_parquet):
     """Point divisions_area at the exonym parquet for the duration of a test."""

tests/test_api.py

@@ -181,6 +181,30 @@ class TestSearchFuzzy:
         assert top["id"] == "manaus"
         assert top["subtype"] == "locality"
 
+    def test_fuzzy_search_puts_the_largest_london_first(
+        self, client, localities_source, urban_centres_source
+    ):
+        # Every candidate matches "London" exactly; population alone orders
+        # them, across sources.
+        resp = client.get(
+            "/search",
+            params={"q": "London", "mode": "fuzzy", "ids_only": "true", "limit": 4},
+        )
+        assert resp.status_code == 200
+        ids = resp.json()["ids"]
+        assert [i["id"] for i in ids] == [
+            "ghsl_5816",
+            "ghsl_2255",
+            "london-oh",
+            "london-ky",
+        ]
+        assert ids[0]["source"] == "urban_centres"
+
+    def test_urban_centre_geometry_by_id(self, client, urban_centres_source):
+        resp = client.get("/geometry/ghsl_5816")
+        assert resp.status_code == 200
+        assert resp.json()["properties"]["name"] == "London"
+
     def test_fuzzy_search_with_sources(self, client):
         try:
             resp = client.get(

tests/test_search.py

@@ -11,6 +11,7 @@ from gazet.search import (
     search_candidates,
     search_divisions_area,
     search_natural_earth,
+    search_urban_centres,
 )
 
 
@@ -321,6 +322,42 @@ class TestLocalities:
         df = get_by_id(con, "manaus", include_geometry=False)
         assert df.iloc[0]["name"] == "Manaus"
 
+    def test_the_more_populous_of_equal_matches_comes_first(
+        self, con, localities_source
+    ):
+        df = search_divisions_area(
+            con, Place(place="London"), limit=2, include_localities=True
+        )
+        assert df["id"].tolist() == ["london-oh", "london-ky"]
+
+
+class TestUrbanCentres:
+    """GHSL urban centres give a polygon to cities Overture has only as a
+    point, such as London; fuzzy search only."""
+
+    def test_the_largest_london_comes_first(self, con, urban_centres_source):
+        df = search_urban_centres(con, Place(place="London"), limit=2)
+        assert df["id"].tolist() == ["ghsl_5816", "ghsl_2255"]
+        assert df.iloc[0]["population"] == 10_408_332
+
+    def test_not_among_default_candidates(self, con, urban_centres_source):
+        ids = [
+            id
+            for df in search_candidates(con, Place(place="London"))
+            for id in df["id"]
+        ]
+        assert "ghsl_5816" not in ids
+
+    def test_an_id_resolves_to_its_source(self, con, urban_centres_source):
+        df = get_by_id(con, "ghsl_5816", include_geometry=False)
+        assert df.iloc[0]["source"] == "urban_centres"
+        assert df.iloc[0]["country"] == "United Kingdom"
+
+    def test_nothing_when_never_built(self, con, monkeypatch):
+        monkeypatch.setattr(search, "URBAN_CENTRES_PATH", "")
+        assert search_urban_centres(con, Place(place="London")).empty
+        assert get_by_id(con, "ghsl_5816", include_geometry=False).empty
+
 
 class TestStoredSearchNames:
     """normalize_geodata stores each record's folded names; a file without