City-Lifeline-Cockpit/tools/city-massing/prepare-district.mjs

// Prepare an auditable district from original rectangular OSM outlines + individually traced roofs.
// No random infill. Exact WGS84/ECEF -> ENU. All heights must be sampled by sample-district.mjs.
import { readFileSync, writeFileSync, mkdirSync } from 'node:fs'
import { resolve } from 'node:path'
import { Cartesian3, Cartographic, Matrix4, Transforms, Math as CM } from 'cesium'
import { PNG } from 'pngjs'
const OUT = resolve('tools/city-massing/district')
mkdirSync(OUT, { recursive: true })
const artifacts = resolve('.scratch/cockpit-shell-panel-layout/verify/artifacts/t35/imagery')
const src = JSON.parse(readFileSync(resolve(artifacts, 'north-source.json')))
const trace = JSON.parse(readFileSync('tools/city-massing/traced-roofs.json'))
const origin = { lon: 110.395, lat: 28.454 }
const mat = Transforms.eastNorthUpToFixedFrame(Cartesian3.fromDegrees(origin.lon, origin.lat))
const inv = Matrix4.inverseTransformation(mat, new Matrix4())
const local = (c) => {
    let a = Matrix4.multiplyByPoint(inv, Cartesian3.fromDegrees(c[0], c[1]), new Cartesian3())
    return [a.x, a.y, a.z]
}
const toGeo = (p) => {
    const c = Cartographic.fromCartesian(Matrix4.multiplyByPoint(mat, new Cartesian3(...p), new Cartesian3()))
    return [CM.toDegrees(c.longitude), CM.toDegrees(c.latitude)]
}
function geoPx([x, y]) {
    const s = 2 ** src.z * 256
    const lon = ((x + src.x0 * 256) / s) * 360 - 180
    const lat = (Math.atan(Math.sinh(Math.PI * (1 - (2 * (y + src.y0 * 256)) / s))) * 180) / Math.PI
    return [lon, lat]
}
const selected = src.features
    .filter((f) => f.ring.length === 5)
    .map((f) => ({
        id: 'osm-' + f.id,
        ring: f.ring.slice(0, 4),
        polygonPx: f.px.slice(0, 4),
        floors: f.id.startsWith('128533681') ? 10 : 5,
        archetype: f.id.startsWith('128533681') ? 'midrise' : 'apartment',
        source: 'OSM original footprint; locally checked against imagery, not a rescaled AABB',
    }))
// Add individually traced roofs only if separate from original controls (traces overlap source excluded below).
for (const f of trace.roofs)
    selected.push({
        ...f,
        ring: f.polygonPx.map(geoPx),
        source: 'individually traced roof on local orthoimagery; floors/facade estimated',
    })
function convexInside(p, poly) {
    let signs = poly.map((a, i) => {
        const b = poly[(i + 1) % poly.length]
        return (b[0] - a[0]) * (p[1] - a[1]) - (b[1] - a[1]) * (p[0] - a[0])
    })
    return signs.every((x) => x >= 0) || signs.every((x) => x <= 0)
}
function overlap(a, b) {
    for (const poly of [a, b])
        for (let i = 0; i < poly.length; i++) {
            const v = poly[i],
                w = poly[(i + 1) % poly.length],
                axis = [w[1] - v[1], v[0] - w[0]]
            const ap = a.map((p) => p[0] * axis[0] + p[1] * axis[1]),
                bp = b.map((p) => p[0] * axis[0] + p[1] * axis[1])
            if (Math.max(...ap) <= Math.min(...bp) || Math.max(...bp) <= Math.min(...ap)) return false
        }
    return true
}
const accepted = [],
    excluded = [],
    samples = []
for (const f of selected) {
    const p = f.ring.map(local)
    const center = p.reduce((a, b) => [a[0] + b[0] / 4, a[1] + b[1] / 4, 0], [0, 0, 0])
    const edges = p.map((a, i) => {
        const b = p[(i + 1) % 4]
        return { len: Math.hypot(b[0] - a[0], b[1] - a[1]), dir: [b[0] - a[0], b[1] - a[1]] }
    })
    const e = edges.reduce((a, b) => (a.len > b.len ? a : b))
    const u = e.dir.map((x) => x / e.len),
        v = [-u[1], u[0]]
    const pro = p.map((a) => [
        (a[0] - center[0]) * u[0] + (a[1] - center[1]) * u[1],
        (a[0] - center[0]) * v[0] + (a[1] - center[1]) * v[1],
    ])
    const width = Math.max(...pro.map((a) => a[0])) - Math.min(...pro.map((a) => a[0]))
    const depth = Math.max(...pro.map((a) => a[1])) - Math.min(...pro.map((a) => a[1]))
    if (width < 5 || depth < 4) {
        excluded.push({ id: f.id, reason: 'too narrow for supported archetype' })
        continue
    }
    // Require a near-rectangle; preserve rotation and original control polygon for audits.
    const rect = [
        [-width / 2, -depth / 2],
        [width / 2, -depth / 2],
        [width / 2, depth / 2],
        [-width / 2, depth / 2],
    ].map(([x, y]) => [center[0] + u[0] * x + v[0] * y, center[1] + u[1] * x + v[1] * y, 0])
    const error = Math.max(...rect.map((q) => Math.min(...p.map((a) => Math.hypot(q[0] - a[0], q[1] - a[1])))))
    if (error > 1.5) {
        excluded.push({ id: f.id, reason: 'irregular outline; needs custom mesh', cornerError: error })
        continue
    }
    if (accepted.some((a) => overlap(a.footprintENU, rect))) {
        excluded.push({ id: f.id, reason: 'overlap with higher-priority source footprint' })
        continue
    }
    const points = []
    const columns = Math.ceil(width / 8),
        rows = Math.ceil(depth / 8)
    for (let y = 0; y <= rows; y++)
        for (let x = 0; x <= columns; x++) {
            const lx = -width / 2 + (width * x) / columns,
                ly = -depth / 2 + (depth * y) / rows
            const q = [center[0] + u[0] * lx + v[0] * ly, center[1] + u[1] * lx + v[1] * ly, 0]
            points.push(toGeo(q))
        }
    const sampleIndices = points.map((g) => {
        samples.push({ lon: g[0], lat: g[1], building: f.id })
        return samples.length - 1
    })
    const centerGeo = toGeo(center)
    accepted.push({
        ...f,
        centerENU: center,
        center: { lon: centerGeo[0], lat: centerGeo[1] },
        footprintENU: rect,
        width,
        depth,
        heading: Math.atan2(u[1], u[0]),
        cornerError: error,
        sampleIndices,
    })
}
const plan = {
    version: 3,
    origin,
    transform: Matrix4.toArray(mat),
    source: trace.description,
    facadesAndFloors: 'design estimates; not survey heights',
    buildings: accepted,
    excluded,
    samples,
}
writeFileSync(resolve(OUT, 'plan.json'), JSON.stringify(plan, null, 2))
const img = PNG.sync.read(readFileSync(resolve(artifacts, 'north-satellite.png')))
for (const f of accepted) {
    const poly = f.polygonPx
    for (let n = 0; n < poly.length; n++) {
        let a = poly[n],
            b = poly[(n + 1) % poly.length]
        const steps = Math.ceil(Math.hypot(b[0] - a[0], b[1] - a[1]))
        for (let i = 0; i <= steps; i++) {
            const x = Math.round(a[0] + ((b[0] - a[0]) * i) / steps),
                y = Math.round(a[1] + ((b[1] - a[1]) * i) / steps)
            if (x < 0 || y < 0 || x >= img.width || y >= img.height) continue
            const k = (y * img.width + x) * 4
            img.data[k] = f.source.startsWith('OSM') ? 255 : 0
            img.data[k + 1] = 220
            img.data[k + 2] = f.source.startsWith('OSM') ? 0 : 255
        }
    }
}
writeFileSync(resolve(OUT, 'footprint-audit.png'), PNG.sync.write(img))
console.log({
    accepted: accepted.length,
    excluded: excluded.length,
    original: accepted.filter((f) => f.id.startsWith('osm-')).length,
    traced: accepted.filter((f) => f.id.startsWith('trace-')).length,
    samplePoints: samples.length,
    out: OUT,
})