City-Lifeline-Cockpit/tools/city-massing/prepare-cad-bim.mjs

// Turn the surveyed CAD footprints into a BIM plan: ENU geometry, tier, height and asset choice.
// Footprints that are close to a rectangle use an approved library asset (windows, balconies, roof
// detail); everything else keeps its real outline as an extruded massing prism, so no building is
// forced into a rectangle it does not have.
import { readFileSync, writeFileSync, mkdirSync } from 'node:fs'
import { resolve } from 'node:path'
import { createHash } from 'node:crypto'
import { Cartesian3, Cartographic, Matrix4, Transforms } from 'cesium'
import * as turf from '@turf/turf'
import { PRIMARY_AREA_BOUNDS, isWithinPrimaryArea } from './primary-area.mjs'

const cityDir = resolve('tools/city-massing/city')
const OUT = resolve(cityDir, 'cad-bim')
mkdirSync(OUT, { recursive: true })

const footprints = JSON.parse(readFileSync(resolve(cityDir, 'cad-buildings.geojson'), 'utf8')).features
const masks = JSON.parse(readFileSync(resolve(cityDir, 'residential-clipped.geojson'))).features
const assetFit = JSON.parse(readFileSync(resolve('tools/city-massing/asset-fit.json'), 'utf8'))
const library = JSON.parse(readFileSync(resolve('tools/city-massing/house-library/library.json'), 'utf8')).archetypes

const origin = { lon: 110.395, lat: 28.454 }
const transform = Transforms.eastNorthUpToFixedFrame(Cartesian3.fromDegrees(origin.lon, origin.lat))
const inverse = Matrix4.inverseTransformation(transform, new Matrix4())
const toLocal = (lon, lat) => Matrix4.multiplyByPoint(inverse, Cartesian3.fromDegrees(lon, lat), new Cartesian3())

// Footprints flatter than this are treated as a rectangle and can host a detailed asset.
const RECT_FILL_MIN = 0.85
const MAX_MASSING_VERTICES = 24

/** Storey estimate from footprint area, in the range this town actually shows. */
function storeysForArea(areaSqm, inMask) {
    if (areaSqm < 40) return 1
    if (areaSqm < 100) return 2
    if (areaSqm < 200) return 3
    if (areaSqm < 400) return 4
    if (areaSqm < 800) return inMask ? 6 : 5
    return inMask ? 8 : 6
}

function orientedRectangle(points) {
    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]
        const length = Math.hypot(dx, dy)
        if (length < 0.05) 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 minU = Math.min(...projected.map((p) => p[0])), maxU = Math.max(...projected.map((p) => p[0]))
        const minV = Math.min(...projected.map((p) => p[1])), maxV = Math.max(...projected.map((p) => p[1]))
        const area = (maxU - minU) * (maxV - minV)
        if (!best || area < best.area)
            best = { area, width: maxU - minU, depth: maxV - minV, heading: Math.atan2(u[1], u[0]) }
    }
    return best
}

function assetFor(width, depth) {
    let best = null
    for (const [family, limits] of Object.entries(assetFit.variants)) {
        const [nominalWidth, nominalDepth] = limits.nominal_m
        for (const [w, d] of [[width, depth], [depth, width]]) {
            const sw = w / nominalWidth, sd = d / nominalDepth
            if (sw < limits.minScale[0] || sw > limits.maxScale[0]) continue
            if (sd < limits.minScale[1] || sd > limits.maxScale[1]) continue
            const penalty = family === 'house' ? 0 : 0.05
            const score = Math.abs(Math.log(sw)) + Math.abs(Math.log(sd)) + penalty
            if (!best || score < best.score) best = { family, score }
        }
    }
    return best ? best.family : null
}

const variants = {
    house: ['house-2f', 'house-4f', 'house'],
    slab: ['slab-4f', 'slab'],
    midrise: ['midrise-8f', 'midrise'],
    'midrise-wide': ['midrise-wide', 'midrise'],
}

const samples = []
const buildings = []
const excluded = []
let assetCount = 0
let massingCount = 0
let outsideMask = 0
let skippedTiny = 0

for (let index = 0; index < footprints.length; index++) {
    const feature = footprints[index]
    if (!isWithinPrimaryArea(feature)) {
        excluded.push({ id: index, reason: 'outside primary area' })
        continue
    }
    const area = turf.area(feature)
    if (area < 4) {
        // Survey noise and structures too small to read as a building at scene scale.
        skippedTiny++
        continue
    }
    const ring = feature.geometry.coordinates[0]
    const local = ring.map(([lon, lat]) => {
        const p = toLocal(lon, lat)
        return [p.x, p.y]
    })
    // Drop the duplicated closing vertex.
    if (local.length > 1 && Math.hypot(local[0][0] - local[local.length - 1][0], local[0][1] - local[local.length - 1][1]) < 0.01)
        local.pop()
    if (local.length < 3) {
        excluded.push({ id: index, reason: 'degenerate outline' })
        continue
    }
    // Simplify heavily for the far field: a 60k-building scene cannot carry full survey outlines.
    const simplified = local.length > MAX_MASSING_VERTICES
        ? turf.simplify(turf.polygon([[...ring]]), { tolerance: 0.00002, highQuality: false }).geometry.coordinates[0]
            .slice(0, -1)
            .map(([lon, lat]) => {
                const p = toLocal(lon, lat)
                return [p.x, p.y]
            })
        : local
    const mask = masks.find((candidate) => turf.booleanIntersects(feature, candidate))
    const inMask = !!mask
    if (!inMask) outsideMask++
    const rect = orientedRectangle(local)
    const fill = rect && rect.area > 0 ? area / rect.area : 0
    const archetype = rect && fill >= RECT_FILL_MIN ? assetFor(rect.width, rect.depth) : null

    const id = `cad-${index}`
    const positions = simplified.map((p) => [Math.round(p[0] * 100) / 100, Math.round(p[1] * 100) / 100])
    // Sample the outline plus the centroid; the seating stage needs real ground under every corner.
    const sampleIndices = [...ring.slice(0, -1), turf.centroid(feature).geometry.coordinates].map(([lon, lat]) => {
        samples.push({ lon, lat, building: id })
        return samples.length - 1
    })

    let kind, width, depth, heading
    if (archetype) {
        kind = 'asset'
        width = rect.width
        depth = rect.depth
        heading = rect.heading
        assetCount++
    } else {
        kind = 'massing'
        width = rect ? rect.width : 6
        depth = rect ? rect.depth : 6
        heading = rect ? rect.heading : 0
        massingCount++
    }
    buildings.push({
        id,
        kind,
        tier: inMask ? 'urban' : 'rural',
        maskId: mask ? mask.properties.ORIG_FID : null,
        footprintENU: positions,
        areaSqm: Math.round(area * 10) / 10,
        rectFillRatio: Math.round(fill * 1000) / 1000,
        width: Math.round(width * 100) / 100,
        depth: Math.round(depth * 100) / 100,
        heading: Math.round(heading * 10000) / 10000,
        archetype,
        preferredAssets: archetype ? variants[archetype] : null,
        floors: storeysForArea(area, inMask),
        heightSource: 'estimated-from-footprint-area',
        sampleIndices,
    })
}

const reasonCounts = {}
for (const entry of excluded) reasonCounts[entry.reason] = (reasonCounts[entry.reason] || 0) + 1

const plan = {
    version: 1,
    source: '1:500 survey CAD layer JMD (居民地), converted through ArcGIS to WGS84',
    policy: 'Real survey outlines are authoritative. Rectangular footprints may host a library asset; all others keep their outline as an extruded massing prism.',
    bounds: [...PRIMARY_AREA_BOUNDS],
    origin,
    transform: Matrix4.toArray(transform),
    buildings,
    samples,
    excluded,
    exclusionReasons: reasonCounts,
    counts: {
        footprints: footprints.length,
        modeled: buildings.length,
        asset: assetCount,
        massing: massingCount,
        urban: buildings.filter((b) => b.tier === 'urban').length,
        rural: buildings.filter((b) => b.tier === 'rural').length,
        skippedTiny,
    },
    heightCaveat: 'Storey counts and heights are estimated from footprint area; the survey carries no usable storey attribute.',
    disclaimer: 'Massing built on surveyed footprints; heights and facades are estimates, not measured.',
}
plan.sha256 = createHash('sha256').update(JSON.stringify(buildings)).digest('hex')
writeFileSync(resolve(OUT, 'plan.json'), JSON.stringify(plan, null, 2))
console.log(JSON.stringify({ ...plan.counts, samples: samples.length, exclusionReasons: reasonCounts, sha256: plan.sha256 }))