# Now Cities methodology

How the bivariate map on [nowcity.org](https://nowcity.org/) is built, what data goes into it, and the modelling choices that shape the picture.

## TL;DR

Take today's global population grid. Subtract the cells that go underwater when every land-ice sheet melts (70 metres of sea level rise). Subtract the cells whose mean annual temperature crosses out of the historical human climate niche under SSP5-8.5. The displaced population is **3.81 billion people**. Push them through a country-level gravity model that combines Geoffrey West's sublinear urban scaling, distance decay, and a WRI Aqueduct future-water-stress penalty. Within each receiving country distribute migrants around existing cities (sharp, point-like, no wide smoothing) plus a frontier share placed as discrete new towns on newly viable cold-out land. Encode the result as a 3 by 3 bivariate raster: existing settlement density on one axis, incoming migration density on the other. Bake to PMTiles and render with MapLibre.

The model is **not a forecast**. It is a counterfactual: *if displaced people distributed by gravity, West-scaling, and water-availability, and then nucleated around existing settlements, what would the picture look like?* Real migration depends on bilateral agreements, conflict, ethnic networks, and a thousand other things this model knows nothing about. The goal here is directional plausibility, not precision.

![Parti diagram synthesizing the whole model: drowned coast, evicted tropics, and poleward resettlement](docs/parti/parti-5-synthesis.svg)

*The whole move at once. Maximum sea level rise drowns the hatched margin, heat and drought evict the tropics (diagonal hatch), and the surviving population plus the new frontier resettle poleward (density tone). The rest of this document explains how each layer is built.*

## 1. The question

> Where do the people go when their home becomes physically uninhabitable?

This map answers it under one specific worst-plausible scenario: max sea level rise (the multi-century all-ice-melted ceiling, 70 m) crossed with max climate warming (SSP5-8.5 at the end of this century). The goal is a defensible high-end planning artifact, not a 2100 forecast.

Other scenarios (lower SLR, lower SSP) produce qualitatively similar but smaller versions of the same picture. The point of pinning to max-max is to give policymakers and urbanists a clear visual answer to the question "where does long-term humanity sit when the climate and the coastlines move to their worst plausible positions?"

![Parti diagram of present-day settlement density](docs/parti/parti-0-today-density.svg)

*Where people live today. Density concentrates in temperate Eurasia, South Asia, and the coasts, fading to empty desert, mountain, and ice. Every later move starts from this baseline.*

## 2. Data sources

All inputs are open and globally consistent. Where projections were needed, the model uses standard CMIP6 forcing under SSP5-8.5.

### 2.1 Population: GHS-POP + GHS-SMOD R2023A

Joint Research Centre Global Human Settlement Layer, 2023A release. Population counts per 1 km cell on the World Mollweide grid, 2020 epoch [^GHS-POP]. Combined with GHS-SMOD R2023A (Settlement Model classification: water, very low density, low density, rural cluster, suburban, low-density urban, dense urban cluster [^GHS-SMOD]).

GHS-POP is the worldwide standard for sub-national population density. For Now Cities it serves three roles: the existing-settlement baseline, the population mass that the gravity model redistributes, and (filtered to SMOD codes 13 and up) the "where people actually live today" mask.

### 2.2 Elevation and bathymetry: GEBCO 2024

GEBCO 2024 global terrain model, sub-ice topography, 15 arc-second resolution [^GEBCO]. Used for the sea-level-rise mask. The sub-ice variant matters because the 70 m ceiling assumes Greenland and Antarctica have fully melted, so the new coastline runs along bedrock, not the present ice surface.

### 2.3 Climate baseline: CHELSA v2.1 bio1

Climatologies at High resolution for the Earth's Land Surface Areas, v2.1, bioclimatic variable 1 (mean annual temperature, MAT). 30 arc-second resolution, baseline period 1981 to 2010 [^CHELSA].

CHELSA's headline win over WorldClim is its better topographic downscaling and superior precipitation. For MAT specifically the two are comparable, but CHELSA's matching CMIP6 future product is what locks us in to this dataset.

### 2.4 Climate future: CHELSA-CMIP6 / GFDL-ESM4

Same CHELSA grid, same bio1 field, but downscaled from GFDL-ESM4 SSP5-8.5 forcing for the 2071 to 2100 horizon [^CHELSA-CMIP6]. GFDL-ESM4 is one of the standard CMIP6 ScenarioMIP runs included in the IPCC AR6 assessment.

We chose a single CMIP6 model (rather than a multi-model ensemble) on purpose. The user-facing scenario is "max-max"; running it through a single warm-biased model produces a clearer worst-case picture than averaging across models would. The price is that we cannot show between-model uncertainty.

### 2.5 Climate niche: Xu et al. 2020

Xu, Kohler, Lenton, Svenning, and Scheffer (2020), *Future of the human climate niche*. PNAS 117(21): 11350-11355 [^Xu2020]. Establishes empirically that the historical distribution of human population is concentrated in a narrow MAT band, with two modes near MAT = 11 to 15 degrees Celsius (the dominant temperate mode) and a small secondary mode near MAT = 25 to 27 degrees Celsius (tropical agricultural belts).

We use Xu's headline finding pragmatically: **MAT > 29 degrees C is the hot-out boundary** (the upper tail beyond which today's humanity essentially does not live). We also extend symmetrically: **MAT < 0 degrees C is the cold-out boundary**, because boreal tundra is below the lower tail of human density in the same data.

A previous iteration used Köppen-Geiger climate classification with "unsuitable" classes (BWh hot desert, EF ice cap, etc.) as the niche mask. That approach was too punishing: 42.9 percent of today's population sat in "bad climate niche" cells under that definition, which violated the model's premise (the niche should describe where people *don't* live). Switching to Xu's MAT thresholds brought today's out-of-niche fraction to 0.9 percent of population, which is the right order of magnitude.

**Planned humid-heat extension.** The current niche cut uses dry MAT alone, with no humidity term. A follow-up version of the model adds a wet-bulb temperature (TWB) cut OR-ed with the MAT cut, drawing on AR6 SYR Fig. 3.2(b) and WG2 Ch.7 (Cissé et al. 2022) which treat combined T+humidity as the governing habitability metric. The target dataset is Copernicus CDS [`sis-extreme-indices-cmip6`](https://cds.climate.copernicus.eu/datasets/sis-extreme-indices-cmip6) annual-max WBT for SSP5-8.5 / 2071-2100; the planned cut is TWB_max > 32 °C as a hard mask, plus a graduated 27-32 °C comfort penalty on existing-city attractor mass per Sherwood & Huber 2010 and Vecellio et al. 2022. The CDS download requires a free ECMWF account and `cdsapi` credentials in `~/.cdsapirc`, so the integration ships in a follow-up version of the model.

### 2.6 Terrain slope: EarthEnv-Topography v1

Pre-computed global slope at 5 km resolution from EarthEnv-Topography v1 [^EarthEnv]. The dataset performs Horn's-algorithm slope on a seamless GMTED2010 mosaic, then aggregates 7.5 arc-sec → 30 arc-sec → 5 arc-min via GRASS and pktools. Output: per-cell **maximum** slope in degrees, single 29.8 MB GeoTIFF, CC-BY 4.0.

We use the max-aggregation variant (`slope_5KMma_GMTEDmd.tif`) so the value at each 5 km cell represents the worst-case slope encountered anywhere inside it. That conservatively excludes cells with steep pockets even if the surrounding area is flat.

Rather than a hard 15° cut, slope feeds the model as a **graduated buildability penalty** on both the existing-city and frontier attractors:

$$
\text{slope\_factor}(x) = \text{clip}\left(\frac{20 - \text{slope}(x)}{12},\ 0,\ 1\right)^{1.5}
$$

| Slope | Factor | Interpretation |
|------:|-------:|----------------|
| ≤ 8° | 1.00 | Hammond plains [^Hammond1964] |
| 12° | 0.60 | Street grids tolerable [^Speck2018] |
| 15° | 0.30 | Hammond hills threshold |
| 18° | 0.10 | Hammond mountains threshold |
| ≥ 20° | 0.00 | Karagulle mountain class [^Karagulle2017] |

The smooth 8°-20° ramp avoids the isobath discontinuities a hard cut produces, and the γ = 1.5 concavity penalises mid-range slopes more aggressively than a linear ramp while keeping the transition continuous. Tunable at the top of `pop-clip-now-cities.py`.

![Parti diagram of steep terrain that resists settlement](docs/parti/parti-3-slopes.svg)

*Steep ground stays unbuilt. Mountain belts above roughly 20 degrees of slope (cross-hatch) resist dense settlement, from the Andes and Rockies to the Himalaya, the Caucasus, and the East African highlands.*

A previous iteration attempted to compute slope ourselves by streaming ~26 k Copernicus GLO-30 30 m COG tiles from AWS, running `gdaldem slope` per tile, and mosaicking via `gdalwarp -r max`. That approach was over-aggressive (3 % global viability vs the ~9 % we should see) and produced visible horizontal scan-line artifacts along the 1° tile boundaries because per-tile bilinear edges didn't line up. Lesson: when a published seamless global slope product exists, use it.

### 2.7 Water stress: WRI Aqueduct 4.0 Future Annual

World Resources Institute Aqueduct 4.0 Water Risk Atlas, Future Annual projection, field **`pes80_ws_x_s`** [^Aqueduct]. Pessimistic emissions (SSP5-8.5 equivalent), 2080 horizon, water-stress score on a 0 to 5 scale at HydroSHEDS Level 6 sub-basins (median basin ~10 000 km²).

Aqueduct's water-stress score is the ratio of total water withdrawals to available renewable supply, with classes:

| Score | Class | Example basins |
|-----:|:------|:---------------|
| 0 | Low | Great Lakes interior, Russian taiga, Patagonia, BC coast, Norway, Brazilian Amazon |
| 1 | Low to medium | US Midwest, eastern Canada, much of Europe, eastern China |
| 2 | Medium to high | US Plains, Mediterranean periphery, Brazilian Cerrado |
| 3 | High | Indian Deccan, South Africa, NSW Australia, Iberia |
| 4 | Extremely high | SW US, North Africa, Punjab, eastern Mediterranean |
| 5 | Catastrophically water-scarce | Arabian Peninsula, Iran, NW China, central Sahara, central Australia |

### 2.7 Future population: UN WPP 2024 (High variant)

United Nations World Population Prospects 2024, total population by country, High variant for 2080 [^WPP2024]. Used to scale today's GHS-POP-derived displaced totals up to 2080 SSP5-equivalent headcounts (§4.5).

The High variant approximates SSP5 at the country aggregate level (global ~14.6 B by 2100 in the High projection, ~12.5 B in published SSP5; close enough as a destination-scaling factor). Multiplier built once via `scripts/wpp-country-multipliers.py`; baseline is the same dataset's Medium variant for 2020 (variants converge for historical years).

### 2.8 Administrative boundaries: Natural Earth

Natural Earth 1:10 million admin-0 country polygons (`ne_10m_admin_0_countries.shp`) [^NaturalEarth]. Public domain. Used to bound the within-country migrant distribution: migrants spread inside a destination country's borders, not across them. The cross-border step is handled by the country-level gravity model.

## 3. Hazard clipping

The base map starts with everyone alive today, then strips two layers of hazard cells.

### 3.1 Sea level rise mask (70 m)

We use 70 m as the **geophysical ice-budget ceiling** for sea level rise — the sum of every land-ice reservoir that could plausibly contribute to SLR if it all eventually melts:

| Source | Sea level equivalent |
|:-------|---------------------:|
| East Antarctic Ice Sheet | ~53 m |
| Greenland Ice Sheet | ~7.4 m |
| West Antarctic Ice Sheet | ~5 m |
| Glaciers and small ice caps | ~0.4 m |
| Antarctic Peninsula | ~0.2 m |
| **Total land ice** | **~66 m** |

The ice reservoirs sum to ~66 m. We plan to **70 m** because the multi-millennial extreme warming required to melt every ice sheet also thermally expands the ocean, adding a few metres of steric rise on top of the ice. Headline figures for a complete land-ice melt already span ~60 m ([NASA](https://sealevel.nasa.gov/understanding-sea-level/global-sea-level/ice-melt/)) to ~70 m ([USGS](https://www.usgs.gov/faqs/how-would-sea-level-change-if-all-glaciers-melted)); 70 m is the high end of the ice budget plus that thermal-expansion term. The 70 m number is the *physics* of the maximum. It is **not** the IPCC's quantified commitment for any specific horizon. AR6 publishes a graded envelope of likely outcomes that sits well below the ice-budget ceiling for any timeline it dares to put numbers on:

| Horizon | SSP5-8.5 likely range | High-impact tail | AR6 confidence |
|:---|---:|:---|:---|
| 2100 (rel. 1995-2014) | 0.63 – 1.01 m | ~2 m if ice-sheet instability triggers | *medium* / *low* tail |
| 2300 | 1.7 – 6.8 m | >15 m low-likelihood storyline | *medium* / *low* + "deep uncertainty" |
| ~4000 CE @ +2 °C | 2 – 6 m | "several metres" from sustained ice loss | *low confidence* |
| Pliocene analogue (~3 Ma, +2.5-4 °C) | n/a | 5 – 25 m higher than today | *medium confidence* |

References: AR6 SYR SPM B.3.1, B.3.3, §3.1.3, Cross-Section Box.2 ([SPM Headline Statements](https://www.ipcc.ch/report/ar6/syr/)).

So 70 m is a planning constant for the question "where on earth do humans have to be out of by the time every land-ice sheet finishes melting." It is not a number any IPCC working group has committed to for 2100 or 2300. We use it because the question of where civilization sits at the *physical limit* of SLR is what Now City asks, and because the answer is robust to its own timing assumptions — a city below 70 m has to move eventually, whether that is 2300 or 8000 CE.

Sub-cell precision: every GEBCO 2024 cell with `elevation <= 70 m above current sea level` becomes ocean. Inland basins below 70 m but disconnected from the ocean (the Caspian, the Salton Sea floor, the Dead Sea) are not flooded by SLR proper. We handle them via a flood-fill from open coast.

The mask is computed from GEBCO 2024 sub-ice elevation: every cell with `elevation <= 70 m above current sea level` becomes ocean. Inland basins that are below 70 m but disconnected from the ocean (the Caspian, the Salton Sea floor, the Dead Sea) are not flooded by SLR proper. We handle them separately with a flood-fill from open coast.

![Parti diagram of the coastline after maximum sea level rise](docs/parti/parti-1-slr.svg)

*Maximum sea level rise redraws the coastline. Seventy metres of melt drowns the hatched margin, taking Hudson Bay, the Gulf and Atlantic coastal plains, Bangladesh, and the West Siberian lowland.*

### 3.2 Climate niche mask (Xu, multi-SSP)

For each climate scenario (today, SSP1-2.6, SSP3-7.0, SSP5-8.5) we classify cells by their MAT against the Xu envelope:

| Today | Future | Class |
|:------|:-------|:------|
| 0 ≤ MAT ≤ 29 °C | 0 ≤ MAT ≤ 29 °C | **stable-in** (viable both periods) |
| 0 ≤ MAT ≤ 29 °C | MAT > 29 °C OR SLR-flooded | **lost** (people need to move) |
| MAT < 0 °C | 0 ≤ MAT ≤ 29 °C AND not flooded | **gained** (newly viable cold-out land) |
| anything else | anything else | **stable-out** (persistently extreme) |

SLR-flooded folds into `lost` regardless of climate class because the person cannot stay either way; the model is indifferent to whether they were evicted by water or by heat.

![Parti diagram of the climate niche shifting poleward](docs/parti/parti-2-niche.svg)

*The climate niche migrates poleward. Heat and drought evict the warm latitudes (diagonal hatch) while newly viable cold-out ground opens up north (solid tone) across Canada, Scandinavia, and Siberia.*

### 3.3 The composite raster

For each `(climate scenario, SLR level)` tuple we bake a composite 2-bit raster encoding `(hazard_class << 2) | populated`:

```
0  stable-in unpop
1  lost unpop
2  gained unpop
3  stable-out unpop
4  stable-in populated
5  lost populated
6  gained populated
7  stable-out populated
```

There are 4 climate scenarios × 7 SLR levels = 28 composites. The Now Cities map locks to the max-max one (`xu_2071_2100_ssp585_slr70m`). The pop-clip detail page lets the user slider through the others.

## 4. Migration model: country-level gravity

The migration step starts from a country-level gravity model. For each origin country `i` (the source of displaced population) and each destination country `j`, the flow is:

$$F_{ij} = P_i^{\text{lost}} \cdot \frac{M_j \cdot \text{bonus}_{ij} / d_{ij}^{\alpha}}{\sum_k M_k \cdot \text{bonus}_{ik} / d_{ik}^{\alpha}}$$

so rows sum to `P_i^lost` (every displaced person ends up somewhere). The non-trivial parts are `M_j`, the bonus terms, and `d_ij`.

### 4.1 Attractor mass (Geoffrey West sublinear scaling)

$$M_j = \left(P_j^{\text{existing}}\right)^{\beta_{\text{infra}}} + \gamma_{\text{new}} \cdot A_j^{\text{gained}}$$

**`P_j^existing`** is the country's current population in cells classified as `stable-in` AND `populated` (GHS-SMOD ≥ 13). This filters out the country's pop that lives in cells that won't be viable in the future (India's lost cells don't become Indian destinations).

**`β_infra = 0.85`** is the Geoffrey West sublinear urban-infrastructure exponent from West (2017), *Scale* [^West2017]. Bigger cities offer more, but with diminishing returns per added capita (housing prices, congestion, land scarcity, visa walls). West's empirical exponent across infrastructure-coupled urban metrics typically falls between 0.83 and 0.87 across the global city-size distribution. Megacities saturate naturally under this rule, which keeps the gravity model from piling everyone into Mumbai or Lagos.

**`γ_new × A_j^gained`** is an additive bonus for countries that gain new habitable land. Virgin land doesn't compound with existing population, it stands on its own. `γ_new = 15 ppl/km²` matches the agricultural-belt density assumption (see §4.5 below).

### 4.2 Distance: pop-weighted centroids

Haversine (great-circle) distance between population-weighted centroids:

- origin centroid `c_i` = pop-weighted mean coordinate of `i`'s lost cells
- destination centroid `c_j` = pop-weighted mean coordinate of `j`'s stable-in-populated cells (falling back to area-weighted gained cells if `j` has no stable-in population)

Geographic centroid (the country's bounding-box centre) is a much worse proxy: it puts Russia's centroid in central Siberia instead of where people actually live (between Moscow and the Urals). Pop-weighted centroids put Russia's attractor near Moscow, where it belongs.

**`α_distance = 1.2`** is a soft distance decay. Voluntary migration follows `α ≈ 2`, but forced (climate) migration historically has weaker distance friction. The Syrian, Sahelian, and Central American diasporas have all extended across continents in spite of distance, particularly when receiving countries opened legal channels. 1.2 sits between voluntary and frictionless.

### 4.3 Cultural and political bonuses

- **`same_country_bonus = 3.0`**. Within-country flows get a 3x pull multiplier. Mumbai-to-Pune stays under one passport, no visa, no language change. Worth tripling.
- **`same_continent_bonus = 1.5`**. Within-continent gets 1.5x.

These multipliers approximate visa friction and diaspora effects without hard-coding bilateral immigration policy. They're a deliberate abstraction: the model knows continent membership, not the USA-Canada free-movement agreement or Schengen.

### 4.4 Zipf-rank penalty on top destinations

After computing `M_j`, the top `zipf_top_n = 5` destinations by mass each have their pull multiplied by `zipf_multiplier = 0.7`.

This is independent of West's scaling. West gives the *baseline* shape; the Zipf penalty is a *political* knob on top. The premise: the global top-5 receiving countries (in this scenario typically Russia, China, India, Canada, USA before water stress) cannot absorb their raw-mass share without political and infrastructural collapse. 0.7 says "the cap reduces each of those by 30 percent and that share gets redistributed to mid-tier destinations."

### 4.5 Scaling displaced population to 2080

The settlement geometry on the map is from GHS-POP 2020, but the model is asking about 2080. Today's GHS-POP times the SSP5-8.5 × 70 m hazard mask gives **3.81 billion people displaced from their current cells**. Two ways to read that number:

1. As the population *of those cells today* that the model has to find homes for.
2. As a stand-in for the 2080 population of those same regions, ignoring growth.

Reading (2) is wrong on the global South in particular, where SSP5 anchors the 21st-century population growth: Nigeria roughly triples (2.6×), Bangladesh 1.7×, India 1.5×, USA 1.5×, while China shrinks (0.77×). We apply per-country growth multipliers from UN World Population Prospects 2024 (High variant, 2080 over 2020) to each country's displaced total before running the gravity model. The empirical settlement *pattern* (where the displaced live today) is preserved; the *headcount* reflects 2080 SSP5 demographics.

After scaling, total displaced under SSP5-8.5 × 70 m × WPP-2024-High-2080 × AR6-hotspot-boost-1.25× comes out to **7.06 B people** (all of which the model routes and places, §5.4), against an AR6-tagged "highly vulnerable" envelope of [3.3 – 3.6 B today](https://www.ipcc.ch/report/ar6/syr/) (SPM A.2.2, *high confidence*) and an AR6 low-elevation coastal exposure projection of [&gt;1 B people by 2050 across all SSPs](https://www.ipcc.ch/report/ar6/syr/) (Cross-Section Box.2). Our number sits well above the AR6 envelopes because it is asking a different question: AR6 counts current compound exposure or near-term coastal flood exposure; ours counts who has to *fundamentally relocate* under the geophysical ice-budget ceiling crossed with end-of-century SSP5 demographics, plus the AR6-recognized compound-risk asymmetry. The model is a counterfactual maximum, not a forecast.

Origin-side scaling is the cleaner mathematical statement than destination-side scaling: 2080 displaced from Bangladesh **is** Bangladesh's 2080 population in lost cells. The two would diverge in countries with mixed origin profiles; we go with origin.

#### 4.5.0 AR6 hotspot displacement boost

After the WPP scaling, an additional 1.25× boost is applied to origin-side lost pop in countries that AR6 explicitly identifies as vulnerability hotspots. The hotspot set follows SYR §2 (line 3762), SPM A.2.5, and WG2 Ch.16: West, Central, and East Africa; South Asia; Central and South America; Caribbean SIDS; South Pacific SIDS; Indian Ocean SIDS (~80 countries total).

The reasoning: the hazard-mask × pop count under-estimates actual displacement in these regions because the model can't see compound risks (conflict, sanitation, water and food insecurity, infrastructure failure) that historically multiply climate displacement. AR6 carries this asymmetry as *high confidence*; we encode it as a 25% boost. Set the multiplier higher (1.5×) to reflect a stronger asymmetry, or to 1.0× to disable.

The combined per-country origin multiplier is `WPP × AR6_boost`:

| Country | WPP 2080/2020 | AR6 boost | Combined |
|:---|---:|---:|---:|
| Bangladesh | 1.72 | 1.25 | 2.15 |
| Nigeria | 2.63 | 1.25 | 3.28 |
| India | 1.48 | 1.25 | 1.85 |
| Indonesia | 1.13 | 1.00 | 1.13 |
| USA | 1.46 | 1.00 | 1.46 |
| Russia | 1.08 | 1.00 | 1.08 |
| China | 0.77 | 1.00 | 0.77 |

### 4.5.1 Capacity sanity check

In parallel to the gravity model we run a global capacity check:

```
need_to_move    = Σ pop_2080 in cells classified as lost
gained_area_km2 = Σ cells classified as gained
gained_capacity = gained_area_km2 × target_density_ppl_km2
deficit         = need_to_move − gained_capacity
```

with `target_density_ppl_km2 = 15`, chosen to match the population density of existing sub-Arctic agricultural belts in northern Canada and Russia (Saskatchewan ~1.7, Manitoba ~2.6, Krasnoyarsk Krai ~1.3, Irkutsk Oblast ~3.0, Russian black-earth region 30 to 60; midpoint of working farms ≈ 15) [^FAO].

This 15 ppl/km² is intentionally **not** the global stable-in average density of ~32. Stable-in includes Western Europe at 200+ ppl/km², and we shouldn't assume the newly-warmed Siberian and Canadian frontier will replicate Western European density within a century.

Under SSP5-8.5 × 70 m × WPP-2024-High-2080 × AR6-hotspot-1.25x:

```
need_to_move      = 7.06 B
gained_area_km2   = 32.33 M km²
gained_capacity   = 485.0 M (at 15 ppl/km²)
deficit           = −6.46 B (short)
```

In words: at agricultural-belt densities the world is **6.5 billion people short** of accommodating max-max displacement on virgin land alone. The gravity model resolves this by routing migrants into existing cities (West sublinear scaling lets cities absorb a lot of extra population) and accepting that some destinations will saturate politically (Zipf penalty) and ecologically (water stress, §4.6).

### 4.6 Water-stress penalty

Per-cell water availability is folded in as a multiplicative penalty:

$$\text{water\_factor}(x) = \left(1 - \text{clip}\left(\frac{\text{ws}(x)}{5}, 0, 0.95\right)\right)^{\gamma_{\text{water}}}$$

with `γ_water = 1.5` (configurable).

| Basin ws | factor | Example regions |
|---------:|-------:|-----------------|
| 0 | 1.00 | Great Lakes, Russian taiga, Patagonia, BC coast, Norway |
| 1 | 0.72 | US Midwest, Eastern Canada, much of Europe |
| 2 | 0.46 | US plains, Mediterranean periphery, Brazilian Cerrado |
| 3 | 0.25 | Indian Deccan, South Africa, NSW Australia, Iberia |
| 4 | 0.09 | SW US, North Africa, Punjab, eastern Mediterranean |
| 5 | 0.01 | Arabian Peninsula, Iran, NW China, central Sahara |

The factor applies **twice**, at different scales:

1. **Country level** in the gravity model: destination mass `M_j` is multiplied by the population-weighted mean `water_factor` over each country's viable populated cells. Saudi Arabia and Iran collapse toward zero; Canada, Russia, Norway, Brazil stay near their baseline West-mass.
2. **Cell level** in the within-country distribution (§5 below): both the existing-city attractor and the frontier attractor are multiplied by the per-cell water factor. Within the USA, Great Lakes cells outweigh Phoenix cells by roughly 30 to 1. Within Australia, the Murray-Darling beats the central desert. Within Russia, the Volga and Ob basins beat the Caspian fringe.

Without the second layer the model would happily concentrate USA's incoming pop in Phoenix and Las Vegas because they're large existing cities. With both layers applied, the picture shifts visibly to the Great Lakes corridor, the Pacific Northwest, and the humid Southeast.

The exponent `γ_water = 1.5` is a moderate prior. Set it lower (`γ = 0.5`) for "water matters but not much"; higher (`γ = 3`) for "climate refugees are explicitly water-seeking." 1.5 keeps basin geography legible without flattening every dry place to invisible.

After all penalties, the top-10 receiving countries under SSP5-8.5 × 70 m × Aqueduct-2080 × WPP-2024-High-2080 are:

| Rank | Country | Incoming (B) |
|----:|:--------|------------:|
| 1 | Russia | 1.58 |
| 2 | China | 0.75 |
| 3 | India | 0.57 |
| 4 | Canada | 0.50 |
| 5 | Nigeria | 0.37 |
| 6 | USA | 0.21 |
| 7 | Indonesia | 0.17 |
| 8 | Pakistan | 0.11 |
| 9 | Norway | 0.10 |
| 10 | Sweden | 0.08 |

Russia leads because it combines low water stress (taiga, Volga, Ob, Lena basins all under ws 1) with massive gained-area (newly viable Siberia under +29 °C). USA is suppressed against its raw West-mass because the Southwest's water stress drags the country mean down. Note that this is **counterfactual** mass redistribution, not a forecast: the model says nothing about whether Russia, China, India, or anyone else would actually open its borders to a billion newcomers.

## 5. Cell-level nucleation (within-country)

The gravity model decides *how many* migrants each country receives. The nucleation step decides *where inside that country* they land.

### 5.1 The point-pattern problem

A naive implementation just smooths the existing-population field with a wide Gaussian kernel (we did this in an earlier draft, σ ≈ 150 km) and distributes incoming migrants proportionally to the smoothed field. The visual result is a series of wide blobs centred on city metros: a smooth wash across the Boston-NY-Philadelphia corridor, a similar wash across Alberta-Saskatchewan-Manitoba.

The user feedback was clear: this doesn't look like settlements. Real human settlement patterns are *pointillistic*: every populated 5 km cell is a discrete dot, separated from its neighbours by empty land. A wash flattens all that into noise.

### 5.2 Two-tier nucleated distribution

The current model splits each country's incoming migrants between two attractors:

1. **Existing-city tier (default share 85 %).** Per-cell weight is `(pop)^1.2` for cells with `pop ≥ 200`, with a tiny σ = 5 km smoothing only to nudge immediate-neighbour cells (so the suburban ring of a city core is reached, but Springfield doesn't bleed into Worcester through empty Berkshire forest). Note that this cell-level exponent is *superlinear*, the opposite of the *sublinear* `β = 0.85` used at the country level (§4.1), and for the opposite reason. A sublinear cell exponent gives small cells disproportionate pull, so incoming migration smears across every populated cell in a continuous settled band; that is what painted European Russia and the central-European plain as a uniform low-density wash rather than a set of distinct cities. A superlinear cell exponent concentrates incoming onto the genuine urban cores (a cell with double the existing population attracts ~2.3 times more migrants, not 1.8), so migration nucleates on discrete cities the way it actually does. Raising the floor to 200 people per cell and tightening the smoothing to 5 km keep each core's halo from merging with its neighbours' into a sheet.
2. **Frontier tier (default share 15 %).** Newly viable (`gained`) cells with `pop < 200` receive the rest, placed as **discrete new towns, not a wash**. An earlier version smoothed the whole gained mask with a small Gaussian, which left the interior of any large gained region at a near-uniform weight and painted a flat sheet of low-density "frontier" across the entire Northwest Territories, northern Manitoba and Quebec, and central Siberia. Real frontier settlement is a lattice of towns with open land between them (look at how the Prairie or Nebraska settled), so the current model builds it in two steps: (a) score each gained cell by closeness to an existing community, a log-compressed version of the present-day viable population surface smoothed at ~60 km (the catchment of a rural service town); (b) take local maxima of that score spaced at least ~80 km apart as new-city centres and drop a ~10 km city footprint at each. New cities seed near existing infrastructure first; the land between hubs gets nothing. It is fine, and realistic, for large stretches of the niche to hold nothing. The deep interior fills last, not first.

The 85/15 split is a prior on the Boston-to-Worcester intuition (real climate migrants pile into existing cities) plus the Krasnoyarsk-to-Yakutsk intuition (some pioneer settlement will open). It's tunable. Raising the frontier share shifts the picture toward more new northern towns; lowering it crams everyone into existing cities. The hub spacing and reach are separate knobs (`HUB_SPACING_CELLS`, `FRONTIER_REACH_SIGMA_CELLS`).

### 5.3 What this looks like

In northern Canada the model concentrates migrants around Edmonton, Calgary, Saskatoon, and the Trans-Canada corridor, plus a scatter of discrete new frontier towns anchored on existing northern communities (Yellowknife, Hay River, Thompson, Churchill) with open land between them. It does *not* paint the territories as a uniform sheet. In Russia migrants concentrate around the Volga axis (Moscow-Nizhny Novgorod-Kazan), the Trans-Siberian corridor (Yekaterinburg-Novosibirsk-Krasnoyarsk), and St. Petersburg, with discrete frontier hubs along the Siberian rivers rather than a wash across the taiga.

![Parti diagram of future settlement density after migration](docs/parti/parti-4-settlement-density.svg)

*Future settlement density. The evicted tropics empty out, and surviving population plus the new frontier concentrate poleward across Canada, northern Europe, and Siberia.*

### 5.4 Conservation: everyone is placed

The distribution conserves population end to end: every person the gravity model assigns to a country is placed inside that country. Two stages used to leak and both are closed. Lost cells that fall outside every Natural Earth polygon (coastal and ocean-edge GHS-POP pixels, ~1-2 % of the total) are now routed to real destinations by mass and distance rather than dropped. And a country that has neither a surviving city nor a frontier hub (it received pull purely from the new-land term on gained land too far from any community) falls back onto its viable cells weighted by buildability, so its incoming still lands. A per-run audit prints the displaced, routed, and placed totals; they now agree to within float rounding (~0.004 %). This conservation is about *placement*, not *capacity*: §4.5.1 is the separate check on whether the gained land can physically hold everyone (it cannot, at 15 ppl/km²).

In the USA the surviving humid Southeast, the Tennessee Valley, the Great Lakes corridor (Cleveland-Detroit-Buffalo-Toronto), and the Pacific Northwest light up with magenta-pop receiving cells. Phoenix, Las Vegas, Albuquerque, and the entire Mojave essentially vanish from the post-warming picture because their water-stress factor crashes the cell-level attractor.

## 6. Stripping abandoned cells

A bias in early renderings: cells where today's population *lives* but where the future climate is `lost` or `stable-out` kept rendering as "existing settlement" in green. The map said "Mumbai is a city" when in fact Mumbai has been entirely evicted by SLR + heat under SSP5-8.5 × 70 m.

The fix: in `pop-clip-now-cities.py`, before quantising the bivariate codes, set `existing_pop = 0` in every `lost` or `stable-out` cell. The post-warming map then renders only the **surviving** existing settlements (the column-0 dark-green ramp) plus incoming migration flows (the orange and magenta cells), with no false signal of "Mumbai still here."

Effects you can verify on the map: Mumbai, Karachi, Dhaka, Bangkok, Manila, Lagos, Kinshasa, Cairo, Khartoum, Baghdad, Riyadh, and most of the Indo-Gangetic plain disappear. The Sahel goes empty. Phoenix, Vegas, and southern Florida are eaten by SLR and lost-niche together.

## 7. Bivariate encoding

The bivariate raster encodes two values per 5 km cell using a 3 × 3 grid. Three bins per axis means 9 codes total, packed in one byte:

```
code = existing_q × 3 + new_q

with existing_q, new_q ∈ {0=empty, 1=town, 2=city}
```

### 7.1 Quantisation thresholds

Per 5 km cell (cell area = 25 km²):

| Bin | Cell threshold (ppl/cell) | Density (ppl/km²) | Reads as |
|----:|--------------------------:|------------------:|:---------|
| 0 | < 200 | < 8 | empty / negligible |
| 1 | 200 – 9 999 | 8 – 400 | town / suburb / rural cluster |
| 2 | ≥ 10 000 | ≥ 400 | dense urban core |

These thresholds are intentionally coarse: the model doesn't usefully resolve finer than this, and the legend stays legible at the size that fits in a sidebar.

### 7.2 Palette: darker = denser, magenta pop at the corner

The 9-colour palette is a classic bivariate-choropleth pattern with a corner-pop trick:

| | new=empty (0) | new=town (1) | new=city (2) |
|:-|:----|:----|:----|
| **existing=city (2)** | dark forest green | dark olive-brown | **MAGENTA POP** |
| **existing=town (1)** | medium leaf green | olive | burnt sienna |
| **existing=empty (0)** | transparent | pale peach | deep rust |

Column 0 alone (no incoming migration) is a pure dark-green settlement ramp. Row 0 alone (no existing settlement) is a pure deep-rust frontier ramp. The off-diagonal blends pick up the natural green + orange = olive / brown mix, going progressively darker.

The far corner (code 8, densest existing × densest new) breaks out of the analogous palette into saturated magenta. This is the standard bivariate-map trick for the "headline" cell: a different hue family ensures the reader's eye lands on the densest-densest cells first. Magenta works because it's the complement of yellow-green at high saturation. We initially tried bright yellow there and it reads too hot, too washed-out at small sizes, and too close to the off-diagonal olives. Magenta is unambiguous.

## 8. Rendering pipeline

### 8.1 GeoTIFF → PMTiles

The full pipeline:

1. Source rasters (GHS-POP, GEBCO, CHELSA) are downsampled and reprojected to a common canonical grid: 5 km cells in EPSG:3857 Web Mercator, extent ±20 037 508 m × ±20 037 508 m, dimensions 8015 × 8015 cells.
2. Hazard composites are built per `(climate, slr)` tuple via `pop-clip-composite.py`. Output: byte-valued GeoTIFFs.
3. Migration classification per `(ssp, slr)` tuple via `pop-clip-migration.py`. Output: byte-valued GeoTIFFs encoding the four shift classes × populated bit.
4. Now Cities bivariate via `pop-clip-now-cities.py`. Output: a single 3 × 3 bivariate raster keyed to the max-max scenario.
5. Each raster goes through `pop-clip-bake-pmtiles.sh`:
   - `gdaldem color-relief` paints the categorical raster into an RGBA GeoTIFF via the per-output `.clr` colour table.
   - `gdal_translate -of MBTILES` generates raster PNG tiles at z0-z6.
   - `gdaladdo` builds the overview pyramid.
   - `pmtiles convert` packages the MBTILES into a single PMTiles archive [^PMTiles].
6. PMTiles archives are uploaded to a public HuggingFace dataset (`datasets/nelsn/now-city-pop-clip`) so the published site can read them via HTTP range requests without paying for tile hosting.

### 8.2 MapLibre + pmtiles protocol

The browser side uses MapLibre GL JS with the `pmtiles` protocol adapter. The Now Cities raster is added as a `raster` source with `tileSize: 512` and rendered with nearest-neighbour resampling so the 5 km bivariate codes stay visibly cell-shaped rather than smoothing into a continuous heatmap.

The basemap is OpenFreeMap's free vector tile basemap (raster output from MapLibre). Country and water features are styled to match the Now City brand palette: deep slate water that visually continues from present-day ocean into post-melt ocean, cream-coloured surviving land.

## 9. What this model is *not*

- **Not a forecast.** It is a counterfactual. Real climate migration depends on bilateral agreements, conflict, ethnic networks, diaspora effects, the political economy of receiving countries, the timing of ice-sheet collapse, and a thousand other factors this model does not encode.
- **Not capacity-tested.** The gravity model places people wherever the attractor pulls. It does not verify that the receiving country has the housing, jobs, water, energy, and political will to absorb that many migrants. The §4.5 capacity check is the global counter-check showing the gap (3.49 B short at 15 ppl/km²) but it does not enforce per-country caps.
- **Not bilateral.** The model does not know that Bangladesh-to-India migration is politically harder than Bangladesh-to-Gulf migration. It only knows distance, mass, and continent. A future refinement could weight by diaspora networks using UN bilateral migration stocks.
- **Not temperature-uncertain.** A single CMIP6 model run (GFDL-ESM4) represents the climate side. The CMIP6 ensemble would give a range, but we deliberately chose the warm-bias single model to keep the worst-case story clear.
- **Not infrastructure-aware.** No road network, no rail network, no airport capacity, no power grid. A serious operational version of this map would weight cell-level attractors by access to transportation and energy infrastructure.

This is just a directional planning artifact for thinking about long-term urbanism in light of climate forcing. It answers the question "where do the people roughly go?" at a coarse enough resolution that the answer is robust to its own assumptions.

## 10. Configuration

All numerical knobs live in `scripts/pop-clip-vectors-config.json` and at the top of `scripts/pop-clip-now-cities.py`. To tune:

1. Edit the JSON or the script constants
2. Re-run `pop-clip-now-cities.py`
3. Re-run `pop-clip-bake-pmtiles.sh`
4. The viewer picks up the new pmtiles automatically (or via cache-busting)

A scholar might want `β_infrastructure = 0.7` (harder city saturation), `γ_new_land = 5` (more pessimistic on virgin land), `α_distance = 1.6` (stricter distance preference), or `γ_water = 3` (water-seeking refugees). The config file makes that one edit plus one rerun.

[^EarthEnv]: Amatulli, G., Domisch, S., Tuanmu, M.-N., Parmentier, B., Ranipeta, A., Malczyk, J., Jetz, W. (2018). *A suite of global, cross-scale topographic variables for environmental and biodiversity modeling*. Scientific Data 5: 180040. [doi:10.1038/sdata.2018.40](https://doi.org/10.1038/sdata.2018.40). Direct download: [https://data.earthenv.org/topography/slope_5KMma_GMTEDmd.tif](https://data.earthenv.org/topography/slope_5KMma_GMTEDmd.tif).

[^Hammond1964]: Hammond, E. H. (1964). *Analysis of properties in land form geography: An application to broad-scale land form mapping*. Annals of the Association of American Geographers.

[^Speck2018]: Speck, J. (2018). *Walkable City Rules: 101 Steps to Making Better Places*. Island Press.

[^Karagulle2017]: Karagulle, D., Frye, C., Sayre, R., Breyer, S., Aniello, P., Vaughan, R., Wright, D. (2017). *Modeling global Hammond landform regions from 250-m elevation data*. Transactions in GIS 21(5): 1040-1060.

[^GHS-POP]: Schiavina, M., Freire, S., Carioli, A., MacManus, K. (2023). *GHS-POP R2023A - GHS population grid multitemporal*. European Commission, Joint Research Centre. [https://human-settlement.emergency.copernicus.eu/](https://human-settlement.emergency.copernicus.eu/)

[^GHS-SMOD]: Schiavina, M., Melchiorri, M., Pesaresi, M. (2023). *GHS- SMOD R2023A - GHS settlement layers, application of the Degree of Urbanisation methodology*. European Commission, JRC. [https://human-settlement.emergency.copernicus.eu/](https://human-settlement.emergency.copernicus.eu/)

[^GEBCO]: GEBCO Compilation Group (2024). *GEBCO 2024 Grid*. [doi:10.5285/1c44ce99-0a0d-5f4f-e063-7086abc0ea0f](https://doi.org/10.5285/1c44ce99-0a0d-5f4f-e063-7086abc0ea0f)

[^CHELSA]: Karger, D. N., Conrad, O., Böhner, J., Kawohl, T., Kreft, H., Soria-Auza, R. W., Zimmermann, N. E., Linder, H. P., Kessler, M. (2017). *Climatologies at high resolution for the Earth's land surface areas*. Scientific Data 4: 170122. [https://chelsa-climate.org/](https://chelsa-climate.org/)

[^CHELSA-CMIP6]: Karger, D. N., Wilson, A. M., Mahony, C., Zimmermann, N. E., Jetz, W. (2020). *Global daily 1km land surface precipitation based on cloud cover-informed downscaling*. Scientific Data 7: 332. CHELSA-CMIP6 product page: [https://chelsa-climate.org/cmip6/](https://chelsa-climate.org/cmip6/)

[^Xu2020]: Xu, C., Kohler, T. A., Lenton, T. M., Svenning, J.-C., Scheffer, M. (2020). *Future of the human climate niche*. Proceedings of the National Academy of Sciences 117(21): 11350-11355. [doi:10.1073/pnas.1910114117](https://doi.org/10.1073/pnas.1910114117)

[^Aqueduct]: Hofste, R., Kuzma, S., Walker, S., Sutanudjaja, E. H., et al. (2019). *Aqueduct 3.0: Updated decision-relevant global water risk indicators*. Technical Note, World Resources Institute. Aqueduct 4.0 release: [https://www.wri.org/data/aqueduct-global-maps-40-data](https://www.wri.org/data/aqueduct-global-maps-40-data). Field: `pes80_ws_x_s` (pessimistic 2080 water stress score).

[^NaturalEarth]: Natural Earth public-domain map data, 1:10m admin-0 countries. [https://www.naturalearthdata.com/](https://www.naturalearthdata.com/)

[^WPP2024]: United Nations, Department of Economic and Social Affairs, Population Division (2024). *World Population Prospects 2024*. Download page: [https://population.un.org/wpp/downloads](https://population.un.org/wpp/downloads). High-variant total population by country, 1950-2100.

[^AR6SYR]: IPCC (2023). *Climate Change 2023: Synthesis Report*. Contribution of Working Groups I, II and III to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change. Core Writing Team, H. Lee and J. Romero (eds.). IPCC, Geneva, Switzerland. [https://www.ipcc.ch/report/ar6/syr/](https://www.ipcc.ch/report/ar6/syr/)

[^West2017]: West, G. (2017). *Scale: The Universal Laws of Life, Growth, and Death in Organisms, Cities, and Companies*. Penguin Press. Sublinear infrastructure scaling exponent typically 0.83 to 0.87.

[^FAO]: FAOSTAT regional population and agricultural land use data. [https://www.fao.org/faostat/](https://www.fao.org/faostat/). Used informally to anchor the 15 ppl/km² agricultural-belt density assumption.

[^PMTiles]: Protomaps PMTiles specification: [https://github.com/protomaps/PMTiles](https://github.com/protomaps/PMTiles). Cloud-optimised single-file vector or raster tile archive with HTTP range requests.

<a id="cities-outlook"></a>

## 11. Cities outlook: a per-city expiry period

The [cities outlook](/cities-outlook.html) page turns the same raster stack into one sentence per large city: when, and by what trigger, the city stops being a place to build for under the max-max scenario. Two clocks run for every city, heat and sea, and both are read at the same three marks: the point at which a tenth, half, and nine-tenths of the city's 2020 population is affected. The expiry period runs from the half mark to the nine-tenths mark of whichever clock reaches the half mark first.

**City set.** Every GHS Urban Centre Database (R2024A) centre with a 2025 population of one million or more, 524 in all. The centre's footprint polygon is the sampling unit; GHS-POP 2020 summed over the footprints is a median 93 percent of the UCDB 2025 figure.

**Heat clock.** For each 1 km cell of the footprint, the local warming per degree of global warming is the cell's SSP5-8.5 change in CHELSA bio1 (2071 to 2100 minus 1981 to 2010, the same GFDL-ESM4 field as §2.4) divided by that model's own global warming between the periods, about 2.8 °C: standard pattern scaling. The global warming level at which the cell's mean annual temperature reaches the 29 °C niche edge from §2.5 is then read off the IPCC AR6 WG1 assessed warming path (0.75 °C for the 1981 to 2010 baseline, 1.09 for 2011 to 2020, then the Table SPM.1 best estimates for 2021 to 2040, 2041 to 2060, and 2081 to 2100; SSP5-8.5 is extended to 2300 with the low end of AR6's 6.6 to 14.1 °C range so later dates stay conservative). Cells are weighted by population and the three marks are the years at which a tenth, half, and nine-tenths of the people live outside the niche, rounded to five years; the SSP3-7.0 path gives a comparison date. Because GFDL-ESM4 is a low-sensitivity model, driving its spatial pattern with the assessed path warms a city roughly 30 percent more by 2090 than the model's own run, which is deliberate: the dates follow the assessed warming, and the raw field values stay in the feed. A check that predicts the SSP3-7.0 field from the SSP5-8.5 pattern lands a quarter of a degree warm on average, the size of error to expect from the scaling. This is temperature only; the wet-bulb extension from §2 would pull the humid coastal cities earlier.

**Sea clock.** For heights from 0 to 70 m, the share of the footprint's population at or below that height on GEBCO 2024 sub-ice elevation at 1 km, counting only cells the world ocean can reach. Reachability runs on a 5 km block-minimum elevation grid seeded from the largest connected component of GHS-SMOD water, so lakes and the Caspian do not seed flooding on their own: Baku at minus 28 m stays dry until the Kuma-Manych sill connects the Caspian at about 26 m, and Death Valley never connects. The half and nine-tenths heights become years with the Max SLR chapter's fastest-plausible schedule, piecewise-linear through 2 m in 2100, 12 m in 2500, 22 m in 3500, and 70 m in 7000. Beside that, each city says how the same height reads in the AR6 Synthesis Report envelope (likely 0.63 to 1.01 m by 2100, about 2 m in the high-impact tail, 1.7 to 6.8 m likely by 2300, 15 m in the 2300 low-likelihood storyline, nothing quantified beyond that). The two are shown side by side and never blended.

**Modifiers.** Aqueduct 2080 water stress over the footprint is a flag (4 and above extremely high, 3 to 4 high), the same prior §4.6 uses; Aqueduct is basin-scale and reads high in some lake-fed basins, so Chicago is flagged. Incoming migrants are the new-density raster from §5 summed over the footprint plus a 5 km ring, and the preference utility (present and future) is the calibrated surface from the attractiveness track.

**Result, 2026-09-09 run.** 267 cities (1.04 billion people today) expire by heat: 13 already sit outside the niche, 35 cross by 2050, 151 between 2050 and 2100, and 68 after 2100 on the extended path. The median heat-expiring city crosses around 2070, and the spread from the 2030s to the 2090s follows how far each city already sits from the edge and how fast its region warms per degree of global warming. 111 (0.43 billion) keep their climate but lose half their people to the sea on the fastest schedule, 16 of them by 2300 and 24 by 2500, the rest between 2500 and 7000. 27 lose a tenth but never half. 119 (0.37 billion) have no expiry date; 73 of those receive incoming migrants worth a quarter or more of their present population, and 57 carry extremely high water stress by 2080. The distribution has two modes, this century by heat and the third millennium and beyond by the sea, with only the lowest deltas and polders in between.

**What it omits.** Defences, subsidence, and storm surge (drowning is a still-water bathtub on an averaged 1 km DEM, which rounds delta cities up by a few metres); humidity; more than one climate model's spatial pattern; and any adaptation term at all. The long-now claim is that adaptation defers a contour but does not move it. Script: `scripts/cities-outlook.py`; feed and schema: `data/cities-outlook/README.md`.

## Reproducibility

All scripts to reproduce this map from source data live in `websites/now_city/scripts/`:

- `pop-clip-pipeline.py` — driver that runs the full pipeline
- `pop-clip-composite.py` — hazard composite per `(climate, slr)` tuple
- `pop-clip-migration.py` — migration classification raster
- `pop-clip-now-cities.py` — bivariate Now Cities raster
- `pop-clip-bake-pmtiles.sh` — bake all GeoTIFFs to PMTiles
- `aqueduct-rasterize.py` — fetch + reproject WRI Aqueduct to 5 km grid

PMTiles archives are mirrored on HuggingFace: [`datasets/nelsn/now-city-pop-clip`](https://huggingface.co/datasets/nelsn/now-city-pop-clip).