diff --git a/server/src/atlas/attractor_matching.rs b/server/src/atlas/attractor_matching.rs
index 8fe664ad3..18a0ba2b2 100644
--- a/server/src/atlas/attractor_matching.rs
+++ b/server/src/atlas/attractor_matching.rs
@@ -18,7 +18,7 @@
 use tracing::{error, warn};
 
 use crate::simulation::generator::{
-    AttractorType, CompatibilityMatrix, GeographicAttractor, SettlementClass,
+    AttractorType, CompatibilityMatrix, GeographicAttractor, SettlementClass, SubBiomeVariant,
 };
 
 // ---------------------------------------------------------------------------
@@ -244,6 +244,8 @@ fn synthetic_attractor(placed: &[CityPlacement], grid_w: u32, grid_h: u32) -> Ge
         position: (row, col),
         attractor_type: AttractorType::PlainCenter,
         strength: 0.5,
+        sub_biome: SubBiomeVariant::TemperateGrassland,
+        terrain_modification_cost: 1.0,
     }
 }
 
@@ -488,6 +490,8 @@ mod tests {
             position: (row, col),
             attractor_type: at,
             strength,
+            sub_biome: SubBiomeVariant::TemperateGrassland,
+            terrain_modification_cost: 1.0,
         }
     }
 
diff --git a/server/src/atlas/features.rs b/server/src/atlas/features.rs
new file mode 100644
index 000000000..58dda1dcf
--- /dev/null
+++ b/server/src/atlas/features.rs
@@ -0,0 +1,615 @@
+//! Geographic feature tag extraction — Layer 1 (D-209).
+//!
+//! After D8 drainage analysis (D-208), this module extracts the 7
+//! `AttractorType` tags from the heightmap + river network. Each attractor has
+//! a pixel position and a normalized `strength` (0.0–1.0) derived from local
+//! terrain quality.
+//!
+//! **D-209 / D-223 reconciliation:** D-209 reads ocean/lake polygons from
+//! `markers.json`, but D-223 reduced markers to a names-only pool — those
+//! polygons no longer exist. Coast and lake cells are therefore derived from
+//! the heightmap itself: ocean = the largest connected below-sea-level water
+//! body; lakes = smaller enclosed below-sea-level bodies.
+//!
+//! **Determinism (D-010 #4):** all collections iterate in sorted/row-major
+//! order; the final attractor list is sorted by `(attractor_type, row, col)`.
+//! No `HashMap`/`HashSet` iteration. `strength` is f32 but is never used as a
+//! sort key.
+
+use std::collections::{HashMap, VecDeque};
+
+use crate::atlas::drainage::DrainageResult;
+use crate::atlas::heightmap::BodyHeightmap;
+use crate::simulation::generator::AttractorType;
+
+/// 8-neighbor offsets (dr, dc). Columns wrap horizontally (equirectangular
+/// globe); rows are bounds-clamped at the poles. Matches `drainage::D8`.
+const NB8: [(i32, i32); 8] = [
+    (-1, 0),
+    (1, 0),
+    (0, 1),
+    (0, -1),
+    (-1, 1),
+    (-1, -1),
+    (1, 1),
+    (1, -1),
+];
+
+/// Minimum spacing (cells) between attractors of an areal type, so coastlines /
+/// valleys / plains yield a sparse, placement-friendly set rather than one
+/// attractor per pixel.
+const MIN_SPACING: i32 = 12;
+
+/// Hard cap on attractors per body (keeps the #955 matching tractable).
+const MAX_ATTRACTORS: usize = 256;
+
+/// A feature before sub-biome classification: position, type, strength.
+/// `layer1` enriches these into `GeographicAttractor`s.
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub struct RawAttractor {
+    pub row: u16,
+    pub col: u16,
+    pub attractor_type: AttractorType,
+    pub strength: f32,
+}
+
+/// Precomputed per-cell terrain fields, shared by feature extraction (D-209)
+/// and sub-biome classification (D-210) so neither recomputes them.
+#[derive(Debug, Clone)]
+pub struct TerrainAnalysis {
+    pub w: usize,
+    pub h: usize,
+    /// `elev < sea_level` (any submerged cell).
+    pub ocean_mask: Vec<bool>,
+    /// Submerged cells not part of the largest water body (enclosed lakes/seas).
+    pub lake_mask: Vec<bool>,
+    /// Chebyshev distance (cells) to the nearest ocean cell or river mouth,
+    /// capped at `WATER_DIST_CAP`. Moisture proxy for habitability/sub-biome.
+    pub water_dist: Vec<u16>,
+    /// Local slope proxy in degrees: `atan(max |Δelev| over 8 neighbors)`.
+    /// Elevation is normalized [0,1]; this is a relative steepness measure.
+    pub slope_deg: Vec<f32>,
+    /// Elevation percentile [0,1] among land cells (ocean cells = 0.0).
+    pub elev_pct: Vec<f32>,
+}
+
+const WATER_DIST_CAP: u16 = 255;
+
+#[inline]
+fn idx(r: usize, c: usize, w: usize) -> usize {
+    r * w + c
+}
+
+#[inline]
+fn wrap_col(c: i32, w: i32) -> usize {
+    c.rem_euclid(w) as usize
+}
+
+impl TerrainAnalysis {
+    /// Compute all shared terrain fields for a body. O(w·h).
+    pub fn analyze(hm: &BodyHeightmap, drainage: &DrainageResult) -> TerrainAnalysis {
+        let w = hm.width as usize;
+        let h = hm.height as usize;
+        let n = w * h;
+        let elev = &hm.data;
+        let sea = hm.sea_level;
+
+        let ocean_mask: Vec<bool> = (0..n).map(|i| elev[i] < sea).collect();
+        let lake_mask = compute_lake_mask(&ocean_mask, w, h);
+        let water_dist = compute_water_dist(&ocean_mask, &drainage.river_network.mouths, w, h);
+        let slope_deg = compute_slope(elev, w, h);
+        let elev_pct = compute_elev_percentile(elev, &ocean_mask, w, h);
+
+        TerrainAnalysis {
+            w,
+            h,
+            ocean_mask,
+            lake_mask,
+            water_dist,
+            slope_deg,
+            elev_pct,
+        }
+    }
+
+    #[inline]
+    pub fn is_ocean(&self, r: usize, c: usize) -> bool {
+        self.ocean_mask[idx(r, c, self.w)]
+    }
+}
+
+/// Largest connected below-sea-level component = ocean; all others = lakes.
+/// Deterministic: BFS seeds scanned row-major; ties broken by lowest cell index.
+fn compute_lake_mask(ocean_mask: &[bool], w: usize, h: usize) -> Vec<bool> {
+    let n = w * h;
+    let mut comp = vec![-1i32; n];
+    let mut comp_sizes: Vec<usize> = Vec::new();
+    let mut next_comp = 0i32;
+
+    for start in 0..n {
+        if !ocean_mask[start] || comp[start] >= 0 {
+            continue;
+        }
+        // Flood fill this component (row-major BFS = deterministic).
+        let mut size = 0usize;
+        let mut q = VecDeque::new();
+        comp[start] = next_comp;
+        q.push_back(start);
+        while let Some(cur) = q.pop_front() {
+            size += 1;
+            let (r, c) = (cur / w, cur % w);
+            for &(dr, dc) in &NB8 {
+                let nr = r as i32 + dr;
+                if nr < 0 || nr >= h as i32 {
+                    continue;
+                }
+                let nc = wrap_col(c as i32 + dc, w as i32);
+                let ni = idx(nr as usize, nc, w);
+                if ocean_mask[ni] && comp[ni] < 0 {
+                    comp[ni] = next_comp;
+                    q.push_back(ni);
+                }
+            }
+        }
+        comp_sizes.push(size);
+        next_comp += 1;
+    }
+
+    if comp_sizes.is_empty() {
+        return vec![false; n]; // no water at all
+    }
+
+    // Largest component (tie → lowest comp id, which is the earliest row-major).
+    let mut ocean_comp = 0i32;
+    let mut best = 0usize;
+    for (cid, &sz) in comp_sizes.iter().enumerate() {
+        if sz > best {
+            best = sz;
+            ocean_comp = cid as i32;
+        }
+    }
+
+    // Lakes = submerged cells in any non-ocean component.
+    (0..n).map(|i| comp[i] >= 0 && comp[i] != ocean_comp).collect()
+}
+
+/// Multi-source BFS Chebyshev distance to nearest ocean cell or river mouth.
+fn compute_water_dist(ocean_mask: &[bool], mouths: &[(u16, u16)], w: usize, h: usize) -> Vec<u16> {
+    let n = w * h;
+    let mut dist = vec![WATER_DIST_CAP; n];
+    let mut q = VecDeque::new();
+    // Seeds in row-major order for determinism.
+    for i in 0..n {
+        if ocean_mask[i] {
+            dist[i] = 0;
+            q.push_back(i);
+        }
+    }
+    for &(mr, mc) in mouths {
+        let i = idx(mr as usize, mc as usize, w);
+        if dist[i] != 0 {
+            dist[i] = 0;
+            q.push_back(i);
+        }
+    }
+    while let Some(cur) = q.pop_front() {
+        let d = dist[cur];
+        if d >= WATER_DIST_CAP {
+            continue;
+        }
+        let (r, c) = (cur / w, cur % w);
+        for &(dr, dc) in &NB8 {
+            let nr = r as i32 + dr;
+            if nr < 0 || nr >= h as i32 {
+                continue;
+            }
+            let nc = wrap_col(c as i32 + dc, w as i32);
+            let ni = idx(nr as usize, nc, w);
+            if dist[ni] > d + 1 {
+                dist[ni] = d + 1;
+                q.push_back(ni);
+            }
+        }
+    }
+    dist
+}
+
+/// Local slope proxy: `atan(max |Δelev| to 8 neighbors)` in degrees.
+fn compute_slope(elev: &[f32], w: usize, h: usize) -> Vec<f32> {
+    let n = w * h;
+    let mut slope = vec![0.0f32; n];
+    for r in 0..h {
+        for c in 0..w {
+            let i = idx(r, c, w);
+            let e = elev[i];
+            let mut max_grad = 0.0f32;
+            for &(dr, dc) in &NB8 {
+                let nr = r as i32 + dr;
+                if nr < 0 || nr >= h as i32 {
+                    continue;
+                }
+                let nc = wrap_col(c as i32 + dc, w as i32);
+                let g = (e - elev[idx(nr as usize, nc, w)]).abs();
+                if g > max_grad {
+                    max_grad = g;
+                }
+            }
+            slope[i] = max_grad.atan().to_degrees();
+        }
+    }
+    slope
+}
+
+/// Elevation percentile [0,1] among land cells; ocean cells get 0.0.
+fn compute_elev_percentile(elev: &[f32], ocean_mask: &[bool], w: usize, h: usize) -> Vec<f32> {
+    let n = w * h;
+    // (scaled_elev, idx) for land cells; integer key for deterministic sort.
+    let mut land: Vec<(i64, usize)> = (0..n)
+        .filter(|&i| !ocean_mask[i])
+        .map(|i| ((elev[i] as f64 * 1_000_000.0) as i64, i))
+        .collect();
+    land.sort_unstable_by(|a, b| a.0.cmp(&b.0).then(a.1.cmp(&b.1)));
+    let mut pct = vec![0.0f32; n];
+    let m = land.len();
+    if m <= 1 {
+        for &(_, i) in &land {
+            pct[i] = 0.5;
+        }
+        return pct;
+    }
+    for (rank, &(_, i)) in land.iter().enumerate() {
+        pct[i] = rank as f32 / (m - 1) as f32;
+    }
+    pct
+}
+
+/// Habitability score [0,1] from elevation band, flatness, and moisture.
+/// Used by `ValleyFloor`/`PlainCenter` strengths and the D-210 classifier.
+pub fn habitability(elev_pct: f32, slope_deg: f32, water_dist: u16) -> f32 {
+    let elev_score = (1.0 - (elev_pct - 0.35).abs() / 0.65).clamp(0.0, 1.0);
+    let flat_score = (1.0 - slope_deg / 15.0).clamp(0.0, 1.0);
+    let water_score = (1.0 - water_dist as f32 / 40.0).clamp(0.0, 1.0);
+    (0.4 * elev_score + 0.3 * flat_score + 0.3 * water_score).clamp(0.0, 1.0)
+}
+
+/// Extract the 7 D-209 attractor tags. Returns raw attractors (no sub-biome),
+/// sorted by `(attractor_type, row, col)` for determinism. Higher-priority
+/// types claim their cells first so a cell is tagged at most once.
+pub fn extract_attractors(
+    hm: &BodyHeightmap,
+    drainage: &DrainageResult,
+    ta: &TerrainAnalysis,
+) -> Vec<RawAttractor> {
+    let w = hm.width as usize;
+    let h = hm.height as usize;
+    let elev = &hm.data;
+    let accum = &drainage.flow_accumulation;
+    let max_accum = drainage.max_accumulation.max(1) as f32;
+
+    let mut claimed = vec![false; w * h];
+    let mut out: Vec<RawAttractor> = Vec::new();
+
+    // O(1) river-cell membership (avoids a binary_search per valley candidate).
+    let mut river_mask = vec![false; w * h];
+    for &(r, c) in &drainage.river_network.river_cells {
+        river_mask[idx(r as usize, c as usize, w)] = true;
+    }
+
+    let claim = |out: &mut Vec<RawAttractor>,
+                 claimed: &mut [bool],
+                 r: usize,
+                 c: usize,
+                 at: AttractorType,
+                 strength: f32| {
+        let i = idx(r, c, w);
+        if claimed[i] {
+            return;
+        }
+        claimed[i] = true;
+        out.push(RawAttractor {
+            row: r as u16,
+            col: c as u16,
+            attractor_type: at,
+            strength: strength.clamp(0.0, 1.0),
+        });
+    };
+
+    // 1. RiverMouth — strength = accum / max_accum.
+    for &(r, c) in &drainage.river_network.mouths {
+        let i = idx(r as usize, c as usize, w);
+        let s = accum[i] as f32 / max_accum;
+        claim(&mut out, &mut claimed, r as usize, c as usize, AttractorType::RiverMouth, s);
+    }
+
+    // 2. RiverCrossing — confluences, strength = accum / max_accum * 0.7.
+    for &(r, c) in &drainage.river_network.confluences {
+        let i = idx(r as usize, c as usize, w);
+        let s = accum[i] as f32 / max_accum * 0.7;
+        claim(&mut out, &mut claimed, r as usize, c as usize, AttractorType::RiverCrossing, s);
+    }
+
+    // 3. CoastalAccess — land within 3 cells of ocean, thinned by spacing.
+    //    strength = 0.6 + coast-density bonus (capped).
+    let mut coastal: Vec<(usize, usize, f32)> = Vec::new();
+    for r in 0..h {
+        for c in 0..w {
+            let i = idx(r, c, w);
+            // water_dist (precomputed) is a cheap pre-filter: only cells within
+            // 3 of water can be coastal, so skip the 49-cell scan for inland.
+            if ta.ocean_mask[i] || claimed[i] || ta.water_dist[i] > 3 {
+                continue;
+            }
+            let near = ocean_cells_within(ta, r, c, 3);
+            if near > 0 {
+                let bonus = (near as f32 / 24.0).min(0.3);
+                coastal.push((r, c, 0.6 + bonus));
+            }
+        }
+    }
+    for (r, c, s) in thin_by_spacing(coastal, &claimed, w) {
+        claim(&mut out, &mut claimed, r, c, AttractorType::CoastalAccess, s);
+    }
+
+    // 4. ValleyFloor — gentle slope, mid elevation, positive habitability.
+    let mut valleys: Vec<(usize, usize, f32)> = Vec::new();
+    for r in 0..h {
+        for c in 0..w {
+            let i = idx(r, c, w);
+            if ta.ocean_mask[i] || claimed[i] {
+                continue;
+            }
+            if river_mask[i] || ta.slope_deg[i] >= 5.0 {
+                continue;
+            }
+            if ta.elev_pct[i] < 0.10 || ta.elev_pct[i] > 0.60 {
+                continue;
+            }
+            let hab = habitability(ta.elev_pct[i], ta.slope_deg[i], ta.water_dist[i]);
+            if hab > 0.0 {
+                valleys.push((r, c, hab));
+            }
+        }
+    }
+    for (r, c, s) in thin_by_spacing(valleys, &claimed, w) {
+        claim(&mut out, &mut claimed, r, c, AttractorType::ValleyFloor, s);
+    }
+
+    // 5. PassEntrance — morphological saddles in higher terrain.
+    //    strength = 1 - elev_pct (lower passes score higher).
+    let mut passes: Vec<(usize, usize, f32)> = Vec::new();
+    for r in 1..h.saturating_sub(1) {
+        for c in 0..w {
+            let i = idx(r, c, w);
+            if ta.ocean_mask[i] || claimed[i] || ta.elev_pct[i] < 0.5 {
+                continue;
+            }
+            if is_saddle(elev, r, c, w, h) {
+                passes.push((r, c, 1.0 - ta.elev_pct[i]));
+            }
+        }
+    }
+    for (r, c, s) in thin_by_spacing(passes, &claimed, w) {
+        claim(&mut out, &mut claimed, r, c, AttractorType::PassEntrance, s);
+    }
+
+    // 6. LakeShore — land adjacent to an enclosed lake. strength = 0.5.
+    let mut shores: Vec<(usize, usize, f32)> = Vec::new();
+    for r in 0..h {
+        for c in 0..w {
+            let i = idx(r, c, w);
+            if ta.ocean_mask[i] || ta.lake_mask[i] || claimed[i] {
+                continue;
+            }
+            if adjacent_to_lake(ta, r, c) {
+                shores.push((r, c, 0.5));
+            }
+        }
+    }
+    for (r, c, s) in thin_by_spacing(shores, &claimed, w) {
+        claim(&mut out, &mut claimed, r, c, AttractorType::LakeShore, s);
+    }
+
+    // 7. PlainCenter — very flat, away from everything. strength = hab * 0.4.
+    let mut plains: Vec<(usize, usize, f32)> = Vec::new();
+    for r in 0..h {
+        for c in 0..w {
+            let i = idx(r, c, w);
+            if ta.ocean_mask[i] || claimed[i] || ta.slope_deg[i] >= 2.0 {
+                continue;
+            }
+            let hab = habitability(ta.elev_pct[i], ta.slope_deg[i], ta.water_dist[i]);
+            plains.push((r, c, hab * 0.4));
+        }
+    }
+    for (r, c, s) in thin_by_spacing(plains, &claimed, w) {
+        claim(&mut out, &mut claimed, r, c, AttractorType::PlainCenter, s);
+    }
+
+    // Cap by keeping the strongest across ALL types (so a coast-heavy body
+    // doesn't starve ValleyFloor/PassEntrance/etc.), then sort the survivors
+    // deterministically by (attractor_type, row, col).
+    if out.len() > MAX_ATTRACTORS {
+        out.sort_by(|a, b| {
+            let sa = (a.strength * 1e6) as i64;
+            let sb = (b.strength * 1e6) as i64;
+            sb.cmp(&sa).then(a.row.cmp(&b.row)).then(a.col.cmp(&b.col))
+        });
+        out.truncate(MAX_ATTRACTORS);
+    }
+    out.sort_by(|a, b| {
+        (a.attractor_type as u8, a.row, a.col).cmp(&(b.attractor_type as u8, b.row, b.col))
+    });
+    out
+}
+
+fn ocean_cells_within(ta: &TerrainAnalysis, r: usize, c: usize, radius: i32) -> usize {
+    let mut count = 0;
+    for dr in -radius..=radius {
+        let nr = r as i32 + dr;
+        if nr < 0 || nr >= ta.h as i32 {
+            continue;
+        }
+        for dc in -radius..=radius {
+            let nc = wrap_col(c as i32 + dc, ta.w as i32);
+            if ta.ocean_mask[idx(nr as usize, nc, ta.w)] {
+                count += 1;
+            }
+        }
+    }
+    count
+}
+
+fn adjacent_to_lake(ta: &TerrainAnalysis, r: usize, c: usize) -> bool {
+    for &(dr, dc) in &NB8 {
+        let nr = r as i32 + dr;
+        if nr < 0 || nr >= ta.h as i32 {
+            continue;
+        }
+        let nc = wrap_col(c as i32 + dc, ta.w as i32);
+        if ta.lake_mask[idx(nr as usize, nc, ta.w)] {
+            return true;
+        }
+    }
+    false
+}
+
+/// Morphological saddle: walking the 8-neighbor ring, the sign of
+/// `(neighbor - cell)` alternates at least 4 times (≥2 higher sectors
+/// separated by ≥2 lower sectors).
+fn is_saddle(elev: &[f32], r: usize, c: usize, w: usize, h: usize) -> bool {
+    // Ring order (clockwise) so transitions are meaningful.
+    const RING: [(i32, i32); 8] = [
+        (-1, 0),
+        (-1, 1),
+        (0, 1),
+        (1, 1),
+        (1, 0),
+        (1, -1),
+        (0, -1),
+        (-1, -1),
+    ];
+    let e = elev[idx(r, c, w)];
+    let mut signs = [0i8; 8];
+    for (k, &(dr, dc)) in RING.iter().enumerate() {
+        let nr = r as i32 + dr;
+        if nr < 0 || nr >= h as i32 {
+            return false; // poles can't be saddles in this scheme
+        }
+        let nc = wrap_col(c as i32 + dc, w as i32);
+        signs[k] = if elev[idx(nr as usize, nc, w)] > e { 1 } else { -1 };
+    }
+    let mut transitions = 0;
+    for k in 0..8 {
+        if signs[k] != signs[(k + 1) % 8] {
+            transitions += 1;
+        }
+    }
+    transitions >= 4
+}
+
+/// Greedy spatial thinning: sort candidates by descending strength (ties by
+/// row, col), keep one per `MIN_SPACING` Chebyshev neighborhood.
+///
+/// Uses a bucket grid (cell size = `MIN_SPACING`) so each candidate only checks
+/// the 3×3 neighboring buckets — O(k) amortized rather than O(k²). The kept
+/// order is fully determined by the sorted candidate iteration, so the internal
+/// `HashMap` (lookup only, never iterated for output) does not affect
+/// determinism.
+fn thin_by_spacing(
+    mut cands: Vec<(usize, usize, f32)>,
+    _claimed: &[bool],
+    _w: usize,
+) -> Vec<(usize, usize, f32)> {
+    // Deterministic order: strength desc, then row, col asc.
+    cands.sort_by(|a, b| {
+        let sa = (a.2 * 1e6) as i64;
+        let sb = (b.2 * 1e6) as i64;
+        sb.cmp(&sa).then(a.0.cmp(&b.0)).then(a.1.cmp(&b.1))
+    });
+    let sp = MIN_SPACING.max(1) as usize;
+    let mut buckets: HashMap<(usize, usize), Vec<(usize, usize)>> = HashMap::new();
+    let mut kept: Vec<(usize, usize, f32)> = Vec::new();
+    for (r, c, s) in cands {
+        let (br, bc) = (r / sp, c / sp);
+        let mut ok = true;
+        'scan: for nbr in br.saturating_sub(1)..=br + 1 {
+            for nbc in bc.saturating_sub(1)..=bc + 1 {
+                if let Some(pts) = buckets.get(&(nbr, nbc)) {
+                    for &(kr, kc) in pts {
+                        let dr = (kr as i32 - r as i32).abs();
+                        let dc = (kc as i32 - c as i32).abs();
+                        if dr.max(dc) < MIN_SPACING {
+                            ok = false;
+                            break 'scan;
+                        }
+                    }
+                }
+            }
+        }
+        if ok {
+            buckets.entry((br, bc)).or_default().push((r, c));
+            kept.push((r, c, s));
+        }
+    }
+    kept
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::atlas::drainage;
+
+    fn slope_grid(w: u32, h: u32) -> Vec<f32> {
+        let n = (w * h) as usize;
+        (0..n)
+            .map(|i| {
+                let r = i / w as usize;
+                let c = i % w as usize;
+                1.0 - (r as f32 / h as f32 * 0.5 + c as f32 / w as f32 * 0.5)
+            })
+            .collect()
+    }
+
+    fn hm(data: Vec<f32>, w: u32, h: u32, sea: f32) -> BodyHeightmap {
+        BodyHeightmap {
+            body_id: "T".into(),
+            width: w,
+            height: h,
+            data,
+            sea_level: sea,
+        }
+    }
+
+    #[test]
+    fn deterministic_extraction() {
+        let h = hm(slope_grid(64, 32), 64, 32, 0.3);
+        let dr = drainage::analyze(&h.data, 64, 32, 0.3);
+        let ta = TerrainAnalysis::analyze(&h, &dr);
+        let a1 = extract_attractors(&h, &dr, &ta);
+        let a2 = extract_attractors(&h, &dr, &ta);
+        assert_eq!(a1, a2, "attractor extraction must be deterministic");
+    }
+
+    #[test]
+    fn attractors_sorted_and_bounded() {
+        let h = hm(slope_grid(128, 64), 128, 64, 0.3);
+        let dr = drainage::analyze(&h.data, 128, 64, 0.3);
+        let ta = TerrainAnalysis::analyze(&h, &dr);
+        let a = extract_attractors(&h, &dr, &ta);
+        assert!(a.len() <= MAX_ATTRACTORS);
+        // Sorted by (type as u8, row, col).
+        for win in a.windows(2) {
+            let ka = (win[0].attractor_type.clone() as u8, win[0].row, win[0].col);
+            let kb = (win[1].attractor_type.clone() as u8, win[1].row, win[1].col);
+            assert!(ka <= kb, "attractors must be sorted");
+        }
+    }
+
+    #[test]
+    fn percentile_in_range() {
+        let h = hm(slope_grid(32, 16), 32, 16, 0.3);
+        let dr = drainage::analyze(&h.data, 32, 16, 0.3);
+        let ta = TerrainAnalysis::analyze(&h, &dr);
+        assert!(ta.elev_pct.iter().all(|&p| (0.0..=1.0).contains(&p)));
+        assert_eq!(ta.slope_deg.len(), 32 * 16);
+    }
+}
diff --git a/server/src/atlas/layer1.rs b/server/src/atlas/layer1.rs
new file mode 100644
index 000000000..7cb090436
--- /dev/null
+++ b/server/src/atlas/layer1.rs
@@ -0,0 +1,186 @@
+//! Layer 1 orchestrator — empty-world topography (#953).
+//!
+//! Runs the full Layer-1 pipeline for one body, in order:
+//!   1. D8 priority-flood drainage (D-208) → river network + basins
+//!   2. shared terrain analysis (ocean/lake masks, water distance, slope,
+//!      elevation percentile) — D-209/D-210 inputs
+//!   3. 7-tag geographic feature extraction (D-209)
+//!   4. sub-biome + terrain_modification_cost classification (D-210)
+//!
+//! Output is the in-memory `Layer1Output`, which maps directly onto
+//! `BodyWorldState` (D-203). Name attachment (D-223) is a separate, cheap step
+//! (`attach_feature_names`) so the compute can be benchmarked in isolation and
+//! names sourced from the DB pool independently.
+//!
+//! **Determinism (D-010 #4):** every stage is deterministic; the same heightmap
+//! yields bit-identical attractors and river networks.
+
+use crate::atlas::body_world_state::{DrainageBasin, RiverNetwork};
+use crate::atlas::drainage::{self, DrainageResult};
+use crate::atlas::features::{self, TerrainAnalysis};
+use crate::atlas::heightmap::BodyHeightmap;
+use crate::atlas::subbiome;
+use crate::simulation::generator::{AttractorType, GeographicAttractor};
+
+/// Full Layer-1 result for one body.
+#[derive(Debug, Clone)]
+pub struct Layer1Output {
+    pub body_id: String,
+    pub river_network: RiverNetwork,
+    pub drainage_basins: Vec<DrainageBasin>,
+    /// Geographic attractors (D-209) with sub-biome + cost (D-210), sorted by
+    /// `(attractor_type, row, col)`.
+    pub attractors: Vec<GeographicAttractor>,
+}
+
+/// Run the Layer-1 topography pipeline for a single body.
+pub fn run_layer1(hm: &BodyHeightmap) -> Layer1Output {
+    let drainage: DrainageResult =
+        drainage::analyze(&hm.data, hm.width, hm.height, hm.sea_level);
+    let ta: TerrainAnalysis = TerrainAnalysis::analyze(hm, &drainage);
+
+    let raw = features::extract_attractors(hm, &drainage, &ta);
+    let attractors: Vec<GeographicAttractor> = raw
+        .iter()
+        .map(|r| {
+            let (sub_biome, terrain_modification_cost) =
+                subbiome::classify(&ta, r.row as usize, r.col as usize);
+            GeographicAttractor {
+                position: (r.row, r.col),
+                attractor_type: r.attractor_type,
+                strength: r.strength,
+                sub_biome,
+                terrain_modification_cost,
+            }
+        })
+        .collect();
+
+    Layer1Output {
+        body_id: hm.body_id.clone(),
+        river_network: drainage.river_network,
+        drainage_basins: drainage.drainage_basins,
+        attractors,
+    }
+}
+
+/// Attach pool names (D-223) to the largest computed rivers and mountains.
+///
+/// Rivers are ranked by mouth strength (a proxy for catchment size) descending;
+/// `RiverMouth` attractors take names from `river_names` in that order. Mountain
+/// names attach to the highest-elevation `Alpine`/`PassEntrance` attractors.
+/// Returns `(river_assignments, mountain_assignments)` as `(position, name)`
+/// pairs; positions that outrun the pool get no name (the pool is finite).
+pub fn attach_feature_names(
+    output: &Layer1Output,
+    river_names: &[String],
+    mountain_names: &[String],
+) -> (Vec<((u16, u16), String)>, Vec<((u16, u16), String)>) {
+    // Rivers: RiverMouth attractors, strongest first (ties by row, col).
+    let mut mouths: Vec<&GeographicAttractor> = output
+        .attractors
+        .iter()
+        .filter(|a| a.attractor_type == AttractorType::RiverMouth)
+        .collect();
+    mouths.sort_by(|a, b| {
+        let sa = (a.strength * 1e6) as i64;
+        let sb = (b.strength * 1e6) as i64;
+        sb.cmp(&sa)
+            .then(a.position.0.cmp(&b.position.0))
+            .then(a.position.1.cmp(&b.position.1))
+    });
+    let rivers = mouths
+        .iter()
+        .zip(river_names.iter())
+        .map(|(a, n)| (a.position, n.clone()))
+        .collect();
+
+    // Mountains: Alpine attractors, strongest first.
+    let mut peaks: Vec<&GeographicAttractor> = output
+        .attractors
+        .iter()
+        .filter(|a| {
+            matches!(
+                a.sub_biome,
+                crate::simulation::generator::SubBiomeVariant::Alpine
+            )
+        })
+        .collect();
+    peaks.sort_by(|a, b| {
+        let sa = (a.strength * 1e6) as i64;
+        let sb = (b.strength * 1e6) as i64;
+        sb.cmp(&sa)
+            .then(a.position.0.cmp(&b.position.0))
+            .then(a.position.1.cmp(&b.position.1))
+    });
+    let mountains = peaks
+        .iter()
+        .zip(mountain_names.iter())
+        .map(|(a, n)| (a.position, n.clone()))
+        .collect();
+
+    (rivers, mountains)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    fn slope_grid(w: u32, h: u32) -> Vec<f32> {
+        let n = (w * h) as usize;
+        (0..n)
+            .map(|i| {
+                let r = i / w as usize;
+                let c = i % w as usize;
+                1.0 - (r as f32 / h as f32 * 0.5 + c as f32 / w as f32 * 0.5)
+            })
+            .collect()
+    }
+
+    fn hm(w: u32, h: u32) -> BodyHeightmap {
+        BodyHeightmap {
+            body_id: "TestBody".into(),
+            width: w,
+            height: h,
+            data: slope_grid(w, h),
+            sea_level: 0.3,
+        }
+    }
+
+    #[test]
+    fn run_layer1_is_deterministic() {
+        let h = hm(128, 64);
+        let o1 = run_layer1(&h);
+        let o2 = run_layer1(&h);
+        assert_eq!(o1.attractors.len(), o2.attractors.len());
+        for (a, b) in o1.attractors.iter().zip(o2.attractors.iter()) {
+            assert_eq!(a.position, b.position);
+            assert_eq!(a.attractor_type, b.attractor_type);
+            assert_eq!(a.strength.to_bits(), b.strength.to_bits());
+            assert_eq!(a.sub_biome, b.sub_biome);
+            assert_eq!(
+                a.terrain_modification_cost.to_bits(),
+                b.terrain_modification_cost.to_bits()
+            );
+        }
+        assert_eq!(
+            o1.river_network.river_cells,
+            o2.river_network.river_cells
+        );
+    }
+
+    #[test]
+    fn produces_attractors_and_costs() {
+        let o = run_layer1(&hm(256, 128));
+        assert!(!o.attractors.is_empty(), "expected some attractors");
+        assert!(o.attractors.iter().all(|a| a.terrain_modification_cost >= 1.0));
+        assert!(o.attractors.iter().all(|a| (0.0..=1.0).contains(&a.strength)));
+    }
+
+    #[test]
+    fn name_attachment_respects_pool_size() {
+        let o = run_layer1(&hm(256, 128));
+        let names = vec!["Aldren".to_string(), "Brook".to_string()];
+        let (rivers, _mtn) = attach_feature_names(&o, &names, &[]);
+        assert!(rivers.len() <= names.len());
+    }
+}
diff --git a/server/src/atlas/mod.rs b/server/src/atlas/mod.rs
index ff0ccbe2b..1e00810a0 100644
--- a/server/src/atlas/mod.rs
+++ b/server/src/atlas/mod.rs
@@ -8,8 +8,11 @@ pub mod block_irregularity;
 pub mod body_world_state;
 pub mod district_mix;
 pub mod drainage;
+pub mod features;
 pub mod gen_queue;
 pub mod heightmap;
+pub mod layer1;
 pub mod rng;
 pub mod skeleton_gen;
+pub mod subbiome;
 pub mod tile_condition;
diff --git a/server/src/atlas/subbiome.rs b/server/src/atlas/subbiome.rs
new file mode 100644
index 000000000..48b001832
--- /dev/null
+++ b/server/src/atlas/subbiome.rs
@@ -0,0 +1,159 @@
+//! Sub-biome variant classification and terrain_modification_cost (D-210).
+//!
+//! Each `GeographicAttractor` (D-209) carries a `SubBiomeVariant` and a
+//! `terrain_modification_cost`. Classification uses four heightmap-derivable
+//! signals (D-210):
+//!   - elevation percentile (of body total) — from `TerrainAnalysis::elev_pct`
+//!   - local slope — `TerrainAnalysis::slope_deg`
+//!   - moisture proxy — distance to nearest river mouth / coast (`water_dist`)
+//!   - temperature proxy — latitude of the equirectangular pixel (`row`)
+//!
+//! `Volcanic` is never emitted here: the Layer-1 inputs carry no volcanic
+//! signal (D-210). It remains in the enum for a future volcanic data source.
+//!
+//! **Determinism:** pure function of integer/float inputs with fixed
+//! thresholds; no RNG, no map iteration. `terrain_modification_cost` is f32
+//! but is never used as a sort key.
+
+use crate::atlas::features::TerrainAnalysis;
+use crate::simulation::generator::SubBiomeVariant;
+
+/// Temperature proxy [0,1] from latitude: 1.0 at the equator (`row == h/2`),
+/// 0.0 at the poles (`row == 0` or `row == h-1`).
+#[inline]
+fn temperature(row: usize, h: usize) -> f32 {
+    if h <= 1 {
+        return 1.0;
+    }
+    let lat = row as f32 / (h - 1) as f32; // 0 = north pole, 1 = south pole
+    1.0 - (lat - 0.5).abs() * 2.0
+}
+
+/// Base infrastructure-build cost per sub-biome (D-210 anchors: grassland 1.0,
+/// coastal lowland 1.4, wetland 3.2, alpine 3.8, volcanic 4.5; the rest
+/// interpolated by buildability).
+fn base_cost(v: SubBiomeVariant) -> f32 {
+    match v {
+        SubBiomeVariant::TemperateGrassland => 1.0,
+        SubBiomeVariant::Savanna => 1.1,
+        SubBiomeVariant::Desert => 1.2,
+        SubBiomeVariant::TemperateForest => 1.3,
+        SubBiomeVariant::CoastalLowland => 1.4,
+        SubBiomeVariant::BorealForest => 1.5,
+        SubBiomeVariant::Tundra => 1.6,
+        SubBiomeVariant::TropicalWet => 2.0,
+        SubBiomeVariant::Wetland => 3.2,
+        SubBiomeVariant::Alpine => 3.8,
+        SubBiomeVariant::Volcanic => 4.5,
+    }
+}
+
+/// Classify the sub-biome and compute `terrain_modification_cost` for the cell
+/// at `(row, col)`. Returns `(variant, cost)`.
+pub fn classify(ta: &TerrainAnalysis, row: usize, col: usize) -> (SubBiomeVariant, f32) {
+    let i = row * ta.w + col;
+    let elev_pct = ta.elev_pct[i];
+    let slope = ta.slope_deg[i];
+    let water_dist = ta.water_dist[i];
+    let temp = temperature(row, ta.h);
+
+    let variant = classify_variant(elev_pct, slope, water_dist, temp);
+
+    // Cost = sub-biome base + a slope surcharge (steeper terrain costs more to
+    // build on), capped so a steep grassland never out-costs flat volcanic.
+    let slope_surcharge = (slope / 12.0).min(1.5);
+    let cost = base_cost(variant) + slope_surcharge;
+
+    (variant, cost)
+}
+
+fn classify_variant(elev_pct: f32, _slope: f32, water_dist: u16, temp: f32) -> SubBiomeVariant {
+    // High elevation dominates → Alpine (mountains, regardless of latitude).
+    if elev_pct > 0.80 {
+        return SubBiomeVariant::Alpine;
+    }
+    // Saturated low ground next to water → Wetland.
+    if water_dist <= 2 && elev_pct < 0.30 {
+        return SubBiomeVariant::Wetland;
+    }
+    // Low ground near a coast → Coastal lowland.
+    if water_dist <= 5 && elev_pct < 0.40 {
+        return SubBiomeVariant::CoastalLowland;
+    }
+    // Cold poleward zones.
+    if temp < 0.20 {
+        return SubBiomeVariant::Tundra;
+    }
+    if temp < 0.40 {
+        return SubBiomeVariant::BorealForest;
+    }
+    // Hot equatorial zones split by moisture.
+    if temp > 0.75 {
+        return if water_dist < 20 {
+            SubBiomeVariant::TropicalWet
+        } else if water_dist < 45 {
+            SubBiomeVariant::Savanna
+        } else {
+            SubBiomeVariant::Desert
+        };
+    }
+    // Temperate mid-latitudes split by moisture.
+    if water_dist > 60 {
+        SubBiomeVariant::Desert
+    } else if water_dist < 25 {
+        SubBiomeVariant::TemperateForest
+    } else {
+        SubBiomeVariant::TemperateGrassland
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::atlas::drainage;
+    use crate::atlas::heightmap::BodyHeightmap;
+
+    fn slope_grid(w: u32, h: u32) -> Vec<f32> {
+        let n = (w * h) as usize;
+        (0..n)
+            .map(|i| {
+                let r = i / w as usize;
+                let c = i % w as usize;
+                1.0 - (r as f32 / h as f32 * 0.5 + c as f32 / w as f32 * 0.5)
+            })
+            .collect()
+    }
+
+    #[test]
+    fn temperature_peaks_at_equator() {
+        assert!((temperature(0, 256) - 0.0).abs() < 0.01);
+        assert!((temperature(255, 256) - 0.0).abs() < 0.01);
+        assert!(temperature(128, 256) > 0.98);
+    }
+
+    #[test]
+    fn alpine_for_high_elevation() {
+        // elev_pct > 0.8 → Alpine regardless of other signals.
+        let (v, cost) = (classify_variant(0.95, 30.0, 100, 0.5), base_cost(SubBiomeVariant::Alpine));
+        assert_eq!(v, SubBiomeVariant::Alpine);
+        assert!(cost > 3.0);
+    }
+
+    #[test]
+    fn classify_is_deterministic_and_bounded() {
+        let h = BodyHeightmap {
+            body_id: "T".into(),
+            width: 64,
+            height: 32,
+            data: slope_grid(64, 32),
+            sea_level: 0.3,
+        };
+        let dr = drainage::analyze(&h.data, 64, 32, 0.3);
+        let ta = TerrainAnalysis::analyze(&h, &dr);
+        let (v1, c1) = classify(&ta, 10, 20);
+        let (v2, c2) = classify(&ta, 10, 20);
+        assert_eq!(v1, v2);
+        assert_eq!(c1.to_bits(), c2.to_bits());
+        assert!(c1 >= 1.0, "cost is at least the grassland baseline");
+    }
+}
diff --git a/server/src/simulation/generator.rs b/server/src/simulation/generator.rs
index 819f0a54b..e756f5b67 100644
--- a/server/src/simulation/generator.rs
+++ b/server/src/simulation/generator.rs
@@ -371,7 +371,7 @@ pub enum TerritorialStatus {
 
 /// The type of terrain feature that attracts settlement placement.
 /// Source: D-195, D-209
-#[derive(Serialize, Deserialize, Clone, Debug, PartialEq, Eq)]
+#[derive(Serialize, Deserialize, Clone, Copy, Debug, PartialEq, Eq)]
 pub enum AttractorType {
     /// Where a river meets sea level or coastline. Historically high-value.
     RiverMouth,
@@ -389,8 +389,29 @@ pub enum AttractorType {
     PlainCenter,
 }
 
+/// Fine-grained terrain classification carried by each `GeographicAttractor`.
+/// Classifies the local terrain more finely than the top-level `SettingType`;
+/// drives ZonePalette modifier selection (D-101) and `terrain_modification_cost`.
+/// Source: D-210
+#[derive(Serialize, Deserialize, Clone, Copy, Debug, PartialEq, Eq)]
+pub enum SubBiomeVariant {
+    TropicalWet,
+    TemperateForest,
+    TemperateGrassland,
+    BorealForest,
+    Tundra,
+    Desert,
+    Savanna,
+    Alpine,
+    Wetland,
+    CoastalLowland,
+    /// No heightmap-derivable signal in the Layer-1 inputs (elevation, slope,
+    /// moisture, latitude); reserved for a future volcanic data source (D-210).
+    Volcanic,
+}
+
 /// A terrain feature at a specific map position that influences city placement scoring.
-/// Source: D-195, D-209
+/// Source: D-195, D-209, D-210
 #[derive(Serialize, Deserialize, Clone, Debug)]
 pub struct GeographicAttractor {
     /// Pixel position in heightmap space [row, col].
@@ -398,6 +419,12 @@ pub struct GeographicAttractor {
     pub attractor_type: AttractorType,
     /// Normalized strength 0.0–1.0. Derived from flow accumulation or habitability score.
     pub strength: f32,
+    /// Fine-grained terrain classification at this position (D-210).
+    pub sub_biome: SubBiomeVariant,
+    /// Infrastructure-build cost multiplier (1.0 = baseline grassland; higher =
+    /// more expensive). Derived from `sub_biome` + local slope. Consumed by the
+    /// attractor-matching pipeline (D-211) to penalize marginal cities. (D-210)
+    pub terrain_modification_cost: f32,
 }
 
 /// Compatibility weights between economic roles and attractor types.