City-Lifeline-Cockpit/tools/city-massing/generate-oldtown-fill.mjs

// Source-backed building plan. A building is only created where an independent footprint source
// says there is one: mapped OSM polygons and imagery-segmented roofs. Road lines are used purely
// as exclusion geometry (never as evidence that a building exists).
import { readFileSync, writeFileSync, mkdirSync, existsSync } from 'node:fs'
import { resolve } from 'node:path'
import { createHash } from 'node:crypto'
import * as turf from '@turf/turf'
import { PRIMARY_AREA_BOUNDS, isWithinPrimaryArea } from './primary-area.mjs'
import { candidateFootprint, isSafeOldTownPlacement } from './oldtown-geometry.mjs'

const dir = resolve('tools/city-massing/city')
const out = resolve(dir, 'oldtown')
mkdirSync(out, { recursive: true })
const masks = JSON.parse(readFileSync(resolve(dir, 'residential-clipped.geojson'))).features
const roads = JSON.parse(readFileSync(resolve(dir, 'roads-clipped.geojson'))).features
const water = JSON.parse(readFileSync(resolve(dir, 'water-clipped.geojson'))).features
const osmBuildings = JSON.parse(readFileSync(resolve(dir, 'osm-buildings-clipped.geojson'))).features
const modeledRegular = JSON.parse(readFileSync(resolve(dir, 'regular-estate-roofs.geojson'))).features
const districtSources = new Set(
    JSON.parse(readFileSync(resolve('tools/city-massing/district/lod-preview/manifest.json')))
        .cells.flatMap((cell) => cell.sourceIds)
)
const imageryPath = resolve(dir, 'imagery-footprints.geojson')
const imagery = existsSync(imageryPath)
    ? JSON.parse(readFileSync(imageryPath)).features
    : []

// Mask reliability, decided by reading the 1:1 audit crops
// (tools/city-massing/city/imagery-footprints-verify.png), not by the segmentation itself.
// Mask 63 is a riverbank sliver dominated by greenhouses and field strips; mask 130 is an
// industry park where bright concrete yards read as roofs. Neither is trustworthy as housing, so
// their imagery footprints are excluded instead of being turned into buildings.
const IMAGERY_TRUSTED_MASKS = new Set([65, 66, 126])

// Carriageway exclusion: a built module may sit beside a street but never on the mapped roadway.
const ROAD_BUFFER_M = 3.5
const roadCorridors = roads.map((road) => {
    try {
        return turf.buffer(road, ROAD_BUFFER_M, { units: 'meters' })
    } catch {
        return null
    }
}).filter(Boolean)

// Families the approved library can carry without distorting facades. Nominal footprints and the
// allowed scale window come from asset-fit.json, the same contract the Blender exporter enforces.
const assetFit = JSON.parse(readFileSync(resolve('tools/city-massing/asset-fit.json'), 'utf8'))
const FAMILIES = Object.entries(assetFit.variants).map(([key, value]) => ({ key, ...value }))
const LATTICE = assetFit.lattice
const VARIANTS = {
    house: ['house-2f', 'house-4f', 'house'],
    slab: ['slab-4f', 'slab'],
    midrise: ['midrise-8f', 'midrise'],
    'midrise-wide': ['midrise-wide', 'midrise'],
}

const roots = []
const excluded = []

const digest = (value) => Number.parseInt(createHash('sha1').update(String(value)).digest('hex').slice(0, 8), 16)

function orientedShape(feature) {
    const ring = feature.geometry.coordinates[0].slice(0, -1)
    const center = ring.reduce((sum, point) => [sum[0] + point[0] / ring.length, sum[1] + point[1] / ring.length], [0, 0])
    const cosLat = Math.cos((center[1] * Math.PI) / 180)
    const metresPerLon = 111320 * cosLat
    const points = ring.map(([lon, lat]) => [lon * metresPerLon, lat * 111320])
    let best = null
    for (let i = 0; i < points.length; i++) {
        const a = points[i], b = points[(i + 1) % points.length]
        const dx = b[0] - a[0], dy = b[1] - a[1], length = Math.hypot(dx, dy)
        if (!length) continue
        const u = [dx / length, dy / length], v = [-u[1], u[0]]
        const projected = points.map((p) => [p[0] * u[0] + p[1] * u[1], p[0] * v[0] + p[1] * v[1]])
        const minX = Math.min(...projected.map((p) => p[0])), maxX = Math.max(...projected.map((p) => p[0]))
        const minY = Math.min(...projected.map((p) => p[1])), maxY = Math.max(...projected.map((p) => p[1]))
        const corners = [[minX, minY], [maxX, minY], [maxX, maxY], [minX, maxY]]
        const error = Math.max(...corners.map((corner) => Math.min(
            ...projected.map((p) => Math.hypot(p[0] - corner[0], p[1] - corner[1]))
        )))
        if (!best || error < best.error)
            best = { center, width: maxX - minX, depth: maxY - minY, heading: Math.atan2(u[1], u[0]), error, u, v, metresPerLon }
    }
    return best
}

function familyFor(width, depth) {
    const ranked = []
    for (const family of FAMILIES) {
        const [nominalAlong, nominalAcross] = family.nominal_m
        for (let rows = 1; rows <= 6; rows++) {
            for (let cols = 1; cols <= 10; cols++) {
                if (rows * cols > LATTICE.maxModulesPerFootprint) continue
                const moduleAlong = (width - 2 * LATTICE.inset_m - LATTICE.gap_m * (cols - 1)) / cols
                const moduleAcross = (depth - 2 * LATTICE.inset_m - LATTICE.gap_m * (rows - 1)) / rows
                if (moduleAlong <= 0 || moduleAcross <= 0) continue
                const alongScale = moduleAlong / nominalAlong
                const acrossScale = moduleAcross / nominalAcross
                if (alongScale < family.minScale[0] || alongScale > family.maxScale[0]) continue
                if (acrossScale < family.minScale[1] || acrossScale > family.maxScale[1]) continue
                // Prefer the family whose natural proportions match the module, then the tiling
                // that uses the fewest, least-stretched modules.
                ranked.push({
                    family, rows, cols, moduleAlong, moduleAcross,
                    alongScale, acrossScale,
                    penalty: Math.abs(Math.log((nominalAlong / nominalAcross) / (moduleAlong / moduleAcross))) * 0.6 +
                        Math.abs(Math.log(alongScale)) + Math.abs(Math.log(acrossScale)) +
                        (rows * cols - 1) * 0.12,
                })
            }
        }
    }
    ranked.sort((a, b) => a.penalty - b.penalty)
    return ranked[0] || null
}

function modulesFor(feature, shape, seed) {
    const width = Math.max(shape.width, shape.depth)
    const depth = Math.min(shape.width, shape.depth)
    const alongIsWidth = shape.width >= shape.depth
    const choice = familyFor(width, depth)
    if (!choice) return null
    const { family, rows, cols } = choice
    const gap = LATTICE.gap_m
    const heading = alongIsWidth ? shape.heading : shape.heading + Math.PI / 2
    // The lattice is laid out inside the source polygon, not just inside its bounding rectangle:
    // the measured outline error is removed from the usable area so modules cannot poke outside.
    const margin = Math.min(2.5, Math.max(LATTICE.inset_m, shape.error + LATTICE.inset_m))
    const usableWidth = width - 2 * margin
    const usableDepth = depth - 2 * margin
    if (usableWidth <= 0 || usableDepth <= 0) return null
    const moduleAlong = (usableWidth - gap * (cols - 1)) / cols
    const moduleAcross = (usableDepth - gap * (rows - 1)) / rows
    if (moduleAlong <= 0 || moduleAcross <= 0) return null
    const alongScale = moduleAlong / family.nominal_m[0]
    const acrossScale = moduleAcross / family.nominal_m[1]
    if (alongScale < family.minScale[0] || alongScale > family.maxScale[0]) return null
    if (acrossScale < family.minScale[1] || acrossScale > family.maxScale[1]) return null
    const modules = []
    for (let row = 0; row < rows; row++) {
        for (let col = 0; col < cols; col++) {
            const alongOffset = -usableWidth / 2 + moduleAlong / 2 + col * (moduleAlong + gap)
            const acrossOffset = -usableDepth / 2 + moduleAcross / 2 + row * (moduleAcross + gap)
            const ux = alongIsWidth ? shape.u[0] : shape.v[0]
            const uy = alongIsWidth ? shape.u[1] : shape.v[1]
            const vx = alongIsWidth ? shape.v[0] : shape.u[0]
            const vy = alongIsWidth ? shape.v[1] : shape.u[1]
            const localX = ux * alongOffset + vx * acrossOffset
            const localY = uy * alongOffset + vy * acrossOffset
            const lon = shape.center[0] + localX / shape.metresPerLon
            const lat = shape.center[1] + localY / 111320
            let footprint = null
            // A module whose corners still escape the source polygon is shrunk before being refused.
            for (const shrink of [0, 0.6, 1.2, 1.8]) {
                const candidate = candidateFootprint({
                    lon, lat, width: moduleAlong - shrink, depth: moduleAcross - shrink, heading,
                })
                if (turf.booleanWithin(candidate, feature)) {
                    footprint = candidate
                    break
                }
            }
            if (!footprint) return null
            modules.push({
                id: `${feature.properties.id}-m${row}-${col}`,
                sourceId: feature.properties.id,
                sourceKind: feature.properties.sourceKind,
                placementSource: feature.properties.sourceKind === 'imagery-segmentation'
                    ? 'imagery-roof-segmentation'
                    : 'osm-building-footprint',
                sourceMask: null,
                sourceEvidence: feature.properties.evidence || 'mapped OSM building polygon',
                roadSource: null,
                roadType: null,
                lon,
                lat,
                width: moduleAlong,
                depth: moduleAcross,
                heading,
                archetype: family.key,
                preferredAssets: VARIANTS[family.key],
                floors: family.key === 'house' ? [2, 3, 4][seed % 3]
                    : family.key === 'slab' ? (seed % 2 ? 4 : 6)
                        : family.key === 'midrise' ? (seed % 2 ? 8 : 10)
                            : 10,
                heightSource: 'estimated',
                sourceFootprintError: shape.error,
                sourcePart: row * cols + col,
                sourceModuleCount: rows * cols,
            })
        }
    }
    return modules
}

function footprintCandidates() {
    const list = []
    for (const feature of osmBuildings) {
        const id = `osm-${feature.properties.osm_id}`
        if (districtSources.has(id)) {
            excluded.push({ id, reason: 'already represented by the audited district sample' })
            continue
        }
        if (feature.geometry.coordinates[0].length - 1 > 8) {
            excluded.push({ id, reason: 'irregular OSM outline; needs custom geometry rather than a rectangle' })
            continue
        }
        list.push({
            ...feature,
            properties: {
                ...feature.properties,
                id,
                sourceKind: 'osm-building',
                evidence: `mapped OSM building polygon ${feature.properties.osm_id}`,
            },
        })
    }
    for (const feature of imagery) {
        if (!IMAGERY_TRUSTED_MASKS.has(feature.properties.maskId)) {
            excluded.push({
                id: feature.properties.id,
                reason: `imagery mask ${feature.properties.maskId} is not a trusted building area (see imagery-footprints-verify.png)`,
            })
            continue
        }
        list.push({
            ...feature,
            properties: {
                ...feature.properties,
                sourceKind: 'imagery-segmentation',
                evidence: `satellite roof segmentation ${feature.properties.id}`,
            },
        })
    }
    return list
}

const candidates = footprintCandidates()
const accepted = []
for (const feature of candidates) {
    const id = feature.properties.id
    if (!isWithinPrimaryArea(feature)) {
        excluded.push({ id, reason: 'footprint outside the configured primary area' })
        continue
    }
    const mask = masks.find((candidate) => turf.booleanWithin(feature, candidate))
    if (!mask) {
        excluded.push({ id, reason: 'footprint is not fully inside a verified residential mask' })
        continue
    }
    if (water.some((polygon) => turf.booleanIntersects(feature, polygon))) {
        excluded.push({ id, reason: 'footprint intersects mapped water' })
        continue
    }
    if (modeledRegular.some((roof) => turf.booleanIntersects(feature, roof))) {
        excluded.push({ id, reason: 'overlaps the separately traced regular estate' })
        continue
    }
    const shape = orientedShape(feature)
    if (!shape || shape.error > 1.75) {
        excluded.push({ id, reason: 'source outline is not rectangular enough for the approved assets' })
        continue
    }
    const modules = modulesFor(feature, shape, digest(id))
    if (!modules) {
        excluded.push({ id, reason: 'approved asset library cannot carry this footprint without distortion' })
        continue
    }
    const placed = []
    let blocked = null
    for (const module of modules) {
        const footprint = candidateFootprint(module)
        if (roadCorridors.some((corridor) => turf.booleanIntersects(footprint, corridor))) {
            blocked = 'module would sit on the mapped roadway'
            break
        }
        // Imagery segments are new geometry, so they must also avoid already mapped OSM buildings.
        const already = feature.properties.sourceKind === 'imagery-segmentation' ? osmBuildings : []
        if (!isSafeOldTownPlacement(footprint, {
            mask,
            roads,
            water,
            existing: already,
            modeled: modeledRegular,
            accepted: accepted.map((entry) => entry.footprint),
        })) {
            blocked = 'module collides with mapped road, water, building or an accepted module'
            break
        }
        placed.push({ module, footprint })
    }
    if (blocked || placed.length !== modules.length) {
        excluded.push({ id, reason: blocked || 'part of the source footprint could not be placed safely' })
        continue
    }
    for (const entry of placed)
        accepted.push({ ...entry.module, sourceMask: `resa-${mask.properties.ORIG_FID}`, footprint: entry.footprint })
    roots.push({
        id,
        sourceKind: feature.properties.sourceKind,
        maskId: mask.properties.ORIG_FID,
        moduleIds: placed.map((entry) => entry.module.id),
        areaSqm: Math.round(turf.area(feature)),
        longSideM: Math.round(Math.max(shape.width, shape.depth) * 10) / 10,
        shortSideM: Math.round(Math.min(shape.width, shape.depth) * 10) / 10,
        headingDeg: Math.round(((shape.heading * 180) / Math.PI) * 10) / 10,
        outlineErrorM: Math.round(shape.error * 100) / 100,
    })
}

const reasonCounts = {}
for (const entry of excluded) reasonCounts[entry.reason] = (reasonCounts[entry.reason] || 0) + 1
const plan = {
    version: 4,
    seed: 20260924,
    scope: 'source-backed urban fill from mapped OSM polygons and satellite roof segmentation',
    bounds: [...PRIMARY_AREA_BOUNDS],
    origin: { lon: 110.395, lat: 28.454 },
    source: 'mapped OSM building polygons + satellite roof segmentation; road/water/mask exclusion',
    policy: 'A building exists only where a footprint source says so. Roads are never building evidence.',
    buildings: accepted.map(({ footprint, ...module }) => module),
    sources: roots,
    excluded,
    exclusionReasons: reasonCounts,
    sourceMasks: masks.map((mask) => mask.properties.ORIG_FID),
    disclaimer:
        'Illustrative massing on source footprints; heights/facades estimated. Missing or irregular source footprints remain unmodeled rather than guessed.',
}
plan.sha256 = createHash('sha256').update(JSON.stringify(plan.buildings)).digest('hex')
writeFileSync(resolve(out, 'plan.json'), JSON.stringify(plan, null, 2))
console.log(JSON.stringify({
    footprintCandidates: candidates.length,
    osmCandidates: osmBuildings.length,
    imageryCandidates: imagery.length,
    sourceBuildings: roots.length,
    modules: plan.buildings.length,
    byKind: {
        osm: roots.filter((r) => r.sourceKind === 'osm-building').length,
        imagery: roots.filter((r) => r.sourceKind === 'imagery-segmentation').length,
    },
    byAsset: Object.fromEntries([...new Set(plan.buildings.map((b) => b.archetype))]
        .map((key) => [key, plan.buildings.filter((b) => b.archetype === key).length])),
    exclusionReasons: reasonCounts,
    sha256: plan.sha256,
}))