City-Lifeline-Cockpit/tools/city-massing/review-large-building-imagery.mjs

// Offline-only TMS evidence. No plan, asset, config, server or browser writes.
import { readFileSync, writeFileSync, existsSync, mkdirSync } from 'node:fs'
import { resolve, join } from 'node:path'
import { PNG } from 'pngjs'
import jpeg from 'jpeg-js'
import * as turf from '@turf/turf'

const out = resolve('.scratch/cockpit-shell-panel-layout/verify/artifacts/t35/large-building-audit')
const root = resolve('tools/city-massing/city/cad-bim')
const read = (p) => JSON.parse(readFileSync(p, 'utf8'))
const audit = read(join(out, 'audit.json'))
const survey = read('tools/city-massing/city/cad-buildings.geojson').features
const seated = new Map(read(join(root, 'seated-plan.json')).buildings.map((b) => [b.id, b]))
const names = ['pilot', 'batch-urban8000-dedup', 'batch-urban8000b', 'batch-urban8000c', 'batch-urban-final', 'batch-urban-tail', 'batch-unfitted-fallback']
const active = new Set(names.flatMap((n) => read(join(root, `lod-${n}-roofed2/manifest.json`)).cells.flatMap((c) => c.buildingIds)))
const gis = process.env.GIS_DIR || 'D:/GIS/GIS'
const z = 19, size = 1792
const px = (lon, lat) => [(lon + 180) / 360 * 2 ** z * 256, (1 - Math.asinh(Math.tan(lat * Math.PI / 180)) / Math.PI) / 2 * 2 ** z * 256]
const windows = [
    ['central-estate', 110.3993, 28.4463, ['cad-20562', 'cad-20556', 'cad-20308', 'cad-20534', 'cad-20301', 'cad-20293', 'cad-20346']],
    ['east-dense-edge', 110.4146, 28.4303, ['cad-3544', 'cad-27407']],
    ['south-central-edge', 110.39955, 28.4260, ['cad-15200']],
    ['southwest-edge', 110.3946, 28.40245, ['cad-14165', 'cad-13681']],
    ['giant-industrial', 110.43065, 28.40835, ['cad-4057', 'cad-4058', 'cad-4060', 'cad-4002', 'cad-4007']],
    ['giant-irregular', 110.4123, 28.4251, ['cad-27730']],
    ['western-estate', 110.3567, 28.4593, ['cad-11179', 'cad-11204', 'cad-11154', 'cad-11167', 'cad-11192', 'cad-11395', 'cad-11386']],
]
// Small bitmap text avoids a font renderer dependency; the CAD IDs are readable in saved PNGs.
const glyph = Object.fromEntries(Object.entries({
    C: ['01110','10001','10000','10000','10000','10001','01110'], A: ['01110','10001','10001','11111','10001','10001','10001'], D: ['11110','10001','10001','10001','10001','10001','11110'], '-': ['00000','00000','00000','11111','00000','00000','00000'],
    0: ['01110','10001','10011','10101','11001','10001','01110'], 1: ['00100','01100','00100','00100','00100','00100','01110'], 2: ['01110','10001','00001','00010','00100','01000','11111'], 3: ['11110','00001','00001','01110','00001','00001','11110'], 4: ['00010','00110','01010','10010','11111','00010','00010'], 5: ['11111','10000','10000','11110','00001','00001','11110'], 6: ['01110','10000','10000','11110','10001','10001','01110'], 7: ['11111','00001','00010','00100','01000','01000','01000'], 8: ['01110','10001','10001','01110','10001','10001','01110'], 9: ['01110','10001','10001','01111','00001','00001','01110'],
}).map(([k, rows]) => [k, rows]))
mkdirSync(out, { recursive: true })
const results = []
for (const [name, lon, lat, labelIds] of windows) {
    const center = px(lon, lat), x0 = Math.floor(center[0] / 256) - 3, y0 = Math.floor(center[1] / 256) - 3
    const image = new PNG({ width: size, height: size }); image.data.fill(255)
    const tiles = [], missingTiles = []
    for (let y = 0; y < 7; y++) for (let x = 0; x < 7; x++) {
        const file = `${gis}/yuanlingPic3/${z}/${x0 + x}/${2 ** z - 1 - y0 - y}.png`
        if (!existsSync(file)) { missingTiles.push(file); continue }
        const bytes = readFileSync(file)
        const tile = bytes[0] === 0xff ? jpeg.decode(bytes) : PNG.sync.read(bytes)
        if (tile.width !== 256 || tile.height !== 256) throw Error(`Unexpected tile size: ${file}`)
        PNG.bitblt(tile, image, 0, 0, 256, 256, x * 256, y * 256)
        tiles.push({ file, payload: bytes[0] === 0xff ? 'jpeg' : 'png' })
    }
    const satellite = `${name}-satellite.png`, overlay = `${name}-outlines.png`
    writeFileSync(join(out, satellite), PNG.sync.write(image))
    const point = (c) => { const p = px(...c); return [p[0] - x0 * 256, p[1] - y0 * 256] }
    const pixel = (x, y, color) => {
        x = Math.round(x); y = Math.round(y)
        if (x < 0 || y < 0 || x >= size || y >= size) return
        const k = (y * size + x) * 4
        image.data[k] = color[0]; image.data[k + 1] = color[1]; image.data[k + 2] = color[2]; image.data[k + 3] = 255
    }
    const line = (a, b, color, width = 1) => {
        const n = Math.max(1, Math.ceil(Math.hypot(b[0] - a[0], b[1] - a[1])))
        for (let i = 0; i <= n; i++) for (let dx = 0; dx < width; dx++) for (let dy = 0; dy < width; dy++) pixel(a[0] + (b[0] - a[0]) * i / n + dx, a[1] + (b[1] - a[1]) * i / n + dy, color)
    }
    const draw = (id, color, width) => {
        const feature = survey[Number(id.slice(4))]
        const rings = feature.geometry.type === 'Polygon' ? feature.geometry.coordinates : feature.geometry.coordinates.flat()
        let intersects = false
        for (const r of rings) {
            const p = r.map(point)
            if (Math.max(...p.map((v) => v[0])) < 0 || Math.min(...p.map((v) => v[0])) >= size || Math.max(...p.map((v) => v[1])) < 0 || Math.min(...p.map((v) => v[1])) >= size) continue
            intersects = true
            for (let i = 1; i < p.length; i++) line(p[i - 1], p[i], color, width)
        }
        return intersects
    }
    let activeOutlines = 0
    for (const id of active) if (draw(id, [40, 255, 110], 1)) activeOutlines++
    const visible = audit.candidates.filter((b) => draw(b.id, [255, 50, 60], labelIds.includes(b.id) ? 2 : 1))
    for (const id of labelIds) {
        const b = visible.find((b) => b.id === id)
        if (!b) continue
        const [cx, cy] = point(b.center), text = id.toUpperCase(), scale = 2
        const x = Math.max(0, Math.min(size - text.length * 12, Math.round(cx))), y = Math.max(0, Math.min(size - 18, Math.round(cy)))
        for (let yy = y - 2; yy < y + 16; yy++) for (let xx = x - 2; xx < x + text.length * 12; xx++) pixel(xx, yy, [15, 15, 15])
        text.split('').forEach((c, i) => glyph[c]?.forEach((row, ry) => [...row].forEach((bit, rx) => {
            if (bit === '1') for (let sy = 0; sy < scale; sy++) for (let sx = 0; sx < scale; sx++) pixel(x + i * 12 + rx * scale + sx, y + ry * scale + sy, [255, 220, 160])
        })))
    }
    writeFileSync(join(out, overlay), PNG.sync.write(image))
    const result = { name, center: [lon, lat], z, x0, y0, width: size, height: size, satellite, overlay, tiles, missingTiles, activeOutlines,
        candidates: visible.map(({ geometry, ...b }) => ({ ...b, labeled: labelIds.includes(b.id), seated: seated.get(b.id) || null })),
        legend: { green: 'active surveyed footprint', red: 'audit candidate; NOT confirmation', label: 'CAD source ID' }, imageryDate: 'unknown; local source has no acquisition timestamp' }
    writeFileSync(join(out, `${name}-evidence.json`), JSON.stringify(result, null, 2))
    results.push({ name, tiles: tiles.length, missingTiles: missingTiles.length, candidates: visible.length, activeOutlines, labeledIds: labelIds.filter((id) => visible.some((b) => b.id === id)) })
    console.log(JSON.stringify(results.at(-1)))
}
writeFileSync(join(out, 'imagery-index.json'), JSON.stringify({ source: `${gis}/yuanlingPic3`, tmsYInverted: true, results }, null, 2))

// Human adjudications below were made by opening the actual PNGs with read_image,
// including unmarked central/western satellite images. Re-review if local imagery changes.
const central = ['cad-20556', 'cad-20308', 'cad-20534', 'cad-20301', 'cad-20293']
const western = ['cad-11179', 'cad-11204', 'cad-11154', 'cad-11167', 'cad-11192', 'cad-11395', 'cad-11386']
const otherLayers = [
    ...read('tools/city-massing/city/regular-estate-roofs.geojson').features,
    ...read('tools/city-massing/district/plan.json').buildings.map((b) => turf.polygon([[...b.ring, b.ring[0]]], { id: b.id })),
]
const intersectionFraction = (a, b) => {
    const i = turf.intersect(turf.featureCollection([turf.feature(a), turf.feature(b)]))
    return i ? turf.area(i) / turf.area(turf.feature(a)) : 0
}
const reviewed = []
function add(id, window, classification, decision, roofMorphology, caveat, exclusionReason = null) {
    const b = audit.candidates.find((b) => b.id === id)
    if (!b) throw Error(`Unknown reviewed ID: ${id}`)
    const s = seated.get(id)
    let candidateOverlap = 0, candidateOverlapId = null
    for (const other of audit.candidates) {
        if (other.id === id || Math.abs(other.center[0] - b.center[0]) > 0.01 || Math.abs(other.center[1] - b.center[1]) > 0.01) continue
        const fraction = intersectionFraction(b.geometry, other.geometry)
        if (fraction > candidateOverlap) { candidateOverlap = fraction; candidateOverlapId = other.id }
    }
    let otherLayerOverlap = 0, otherLayerOverlapId = null
    for (const f of otherLayers) {
        const fraction = intersectionFraction(b.geometry, f.geometry)
        if (fraction > otherLayerOverlap) { otherLayerOverlap = fraction; otherLayerOverlapId = f.properties.id }
    }
    reviewed.push({ id, locationWGS84: b.center, areaSqm: b.areaSqm, dimensionsM: [b.width, b.depth], sourceTier: b.tier,
        classification, decision, roofMorphology, imageryMatched: classification !== 'erroneous-compound-outline', imageryReviewed: true,
        evidence: [`${window}-satellite.png`, `${window}-outlines.png`, `${window}-evidence.json`], auditStatus: b.status,
        overlapActiveFraction: b.overlapFraction, candidateOverlapFraction: +candidateOverlap.toFixed(6), candidateOverlapId,
        otherNeighborhoodLayerOverlapFraction: +otherLayerOverlap.toFixed(6), otherNeighborhoodLayerOverlapId: otherLayerOverlapId,
        seating: s ? { alreadySeated: true, construction: s.construction, reliefM: s.terrainRelief, terrainMinM: s.terrainMin,
            terrainMaxM: s.terrainMax, positionENU: s.positionENU, estimatedHeightM: s.heightM, heightSource: s.heightSource } : { alreadySeated: false, reliefM: b.terrainRelief },
        caveat, exclusionReason })
}
for (const id of central) add(id, 'central-estate', 'residential-slab', 'recommend-pilot-flat',
    'Long gray roof bar in a dense riverside residential group; visible ridge/eave line, repeated small roof details and lane separations. The CAD edge follows the roof bar; cad-20308 has genuine stepped/notched outline, not a bounding rectangle.',
    'Rural source tier is a mask/selection classification, not evidence of rural land use. Roof-to-ground parallax and undated imagery remain; floors/facades are estimates. Preserve current seated footprint and gray roof, not a generic blue warehouse roof.')
add('cad-20562', 'central-estate', 'residential-slab', 'recommend-pilot-bounded-plinth',
    'Third long gray parallel residential roof south of the river road; outline follows a separate roof and stops before its active green neighbor.',
    'Already seated with 3.18m DEM relief, only 0.02m below the 3.2m relief cap. Keep the existing bounded plinth; do not widen/raise it. Prefer the five flat candidates for a stricter first export.')
for (const id of western) add(id, 'western-estate', 'residential-slab', 'recommend-second-pilot-flat',
    'Repeated narrow gray pitched/ridged roof bars with light eave edges in a planned dense estate, arranged in parallel rows with real access lanes and a basketball court. Individual CAD outlines match individual roof bars and visible small edge notches.',
    'This is a substantial dense residential estate despite rural source tier; do not restore surrounding scattered houses wholesale. Repetition supports shared detail treatment, not invented identical measured heights. Small source roof features may be blurred.')
for (const id of ['cad-4057', 'cad-4058', 'cad-4060']) add(id, 'giant-industrial', 'industrial-roof', 'exclude-residential-pilot',
    'Huge blue metal-sheet hall roofs with long central ridge and repeated transverse seams; surrounding loading/access strips. CAD outlines match the real industrial halls, not bogus giant GIS rectangles. cad-4060 is only partially inside this mosaic.',
    'Real missing industrial coverage, not missing residential slabs. Not already seated; DEM relief exceeds cap. Use a later purpose-built low-rise industrial model with independent terrain solution, not tall apartment floors.',
    'Wrong land-use/morphology for housing pilot and rejected DEM seating; do not relax the global plinth gate.')
add('cad-27730', 'giant-irregular', 'erroneous-compound-outline', 'reject-wholesale-extrusion',
    'A giant irregular red boundary spans bare brown construction yard, separate bright roofs, courtyards and gaps along the river. No continuous roof matches its complete polygon.',
    'Small true roofs inside this compound may merit later independent IDs; this parent outline is not one building. The 39.19m relief is a warning consistent with an invalid building-level envelope.',
    'Compound/erroneous GIS envelope rather than one roof; would fabricate a giant roof over open ground.')
add('cad-20346', 'central-estate', 'duplicate-survey-outline', 'exclude-pending-dedup',
    'The red square roof east of the residential bars has nested red survey outlines, not two separate structures. cad-20346 is entirely contained by cad-20347.',
    'The audit only compared candidates with active footprints, so zero active overlap does not establish candidate uniqueness. cad-20347 also overlaps its smaller outline by 76.85% of its own area; choose exactly one after roof/structure review.',
    '100% contained in cad-20347; adding both would duplicate one roof. Neither belongs in the selected residential bar pilot.')
add('cad-14165', 'southwest-edge', 'industrial-or-service-roof', 'exclude-residential-pilot',
    'Broad green shed-like rectangular roof in a materials yard with nearby blue roofs, a narrow river and sparse red-roof houses; no repeated apartment row morphology.',
    'Already seated at 3.18m relief does not prove residential use or justify its estimated 19.2m height. Sparse vicinity is less compelling than the confirmed dense estates.',
    'Not compelling residential infill; broad service/industrial roof and uncertain height.')
add('cad-3544', 'east-dense-edge', 'industrial-or-service-roof', 'exclude-residential-pilot',
    'Single plain blue-gray elongated roof on east riverbank amid disturbed/material-storage ground; sports/campus buildings sit across the river, not in the same housing row.',
    'Roof and survey agree and existing relief is only 1.78m, but safe DEM seating is not evidence of apartment use; estimated 19.2m height is unsupported.',
    'Service/industrial-looking roof outside dense apartment rows; retain for separate use/height review.')
const primaryIds = central
const extendedIds = [...central, 'cad-20562']
const secondaryIds = western
for (const b of reviewed.filter((b) => b.decision.startsWith('recommend'))) {
    if (!b.seating.alreadySeated || b.overlapActiveFraction > 0.01 || b.candidateOverlapFraction > 0.05 || b.otherNeighborhoodLayerOverlapFraction > 0.05) throw Error(`Pilot uniqueness/seating failed: ${b.id}`)
}
const report = {
    version: 1, sourceAuditPlanSha256: audit.sourcePlanSha256, activeCoverageDefinition: 'Seven roofed2 CAD manifests, matching audit-large-building-coverage.mjs; district and regular-estate survey/trace layers additionally checked for pilot overlap.',
    method: 'Offline local TMS z19, 7x7 tiles per view, 1792px cap; PNG or JPEG payload decoded from .png files. North-up XYZ rows read with TMS inverted Y. Actual saved images visually reviewed with read_image; no generated/invented satellite imagery.',
    imagery: { directory: `${gis}/yuanlingPic3`, zoom: z, totalTiles: results.reduce((n, r) => n + r.tiles, 0), missingTiles: results.reduce((n, r) => n + r.missingTiles, 0), acquisitionDate: 'unknown', groundResolutionMApprox: 0.263, windows: results },
    findings: 'Confirmed dense residential omissions exist in safely seated rural-tier IDs. Top-area candidates are not interchangeable: real industrial halls and erroneous compound outlines both occur. Zero active overlap alone misses duplicates among nonactive candidates.',
    pilot: { primaryIds, extendedIds, secondaryIds, primaryRationale: 'Five real roof-matched missing residential buildings, all already seated flat, all zero active/candidate/other-neighborhood overlap. Restore source outlines without changing global seating gates.',
        extendedRationale: 'cad-20562 adds a sixth confirmed housing slab but retains its existing near-cap bounded plinth.', secondaryRationale: 'Seven already-seated flat dense-estate slabs with matching gray ridged roofs; a second small batch, not all rural candidates.' },
    candidates: reviewed,
    limits: ['Missing is relative to checked source/model footprint layers, not browser runtime confirmation.', 'Imagery date, measured heights, precise floors, facade colors and windows remain unknown; keep all height/facade metadata explicitly estimated.', 'Do not use unmasked satellite pixels as facade textures. Roof detail must preserve ridge/eave orientation; nearby shadows/trees/roads must not be projected onto roofs.', 'No source plan, assets or config modified; recommendations only. Build/dev/server/browser debug intentionally skipped under AGENTS.md.'],
}
writeFileSync(join(out, 'review.json'), JSON.stringify(report, null, 2))
const artifact = '.scratch/cockpit-shell-panel-layout/verify/artifacts/t35/large-building-audit/'
const link = (file, label = file) => `[${label}](<${artifact}${file}>)`
const rows = reviewed.map((b) => `| ${b.id} | ${b.locationWGS84.map((v) => v.toFixed(7)).join(', ')} | ${b.areaSqm} | ${b.seating.reliefM} / ${b.seating.alreadySeated ? b.seating.construction : 'rejected'} | ${b.decision} | ${link(b.evidence[1], 'image')} |`)
const details = reviewed.map((b) => `### ${b.id}\n- Location (lon, lat): ${b.locationWGS84.join(', ')}; ${b.areaSqm} m².\n- Roof: ${b.roofMorphology}\n- Decision: **${b.decision}**. ${b.exclusionReason || 'Confirmed roof match and already seated; recommended bounded export.'}\n- Evidence: ${b.evidence.map((file) => link(file)).join(' · ')}\n- Caveat: ${b.caveat}\n- Overlap: active ${(b.overlapActiveFraction * 100).toFixed(2)}%; other candidate ${(b.candidateOverlapFraction * 100).toFixed(2)}% (${b.candidateOverlapId || 'none'}); district/regular-estate ${(b.otherNeighborhoodLayerOverlapFraction * 100).toFixed(2)}%.\n- Seating: ${b.seating.alreadySeated ? `${b.seating.construction}, relief ${b.seating.reliefM}m; height ${b.seating.estimatedHeightM}m is estimated, not surveyed` : `not seated; DEM relief ${b.seating.reliefM}m`}.`)
writeFileSync(join(out, 'review.md'), `# Issue 35 — evidence-backed large-building imagery review\n\n## Recommended small pilot\n\n**Primary flat-only five:** ${primaryIds.join(', ')}. These are actual missing residential roof bars in a dense riverside block; all already seated, all unique against checked active CAD, nonactive candidates, district and regular-estate layers. **Optional sixth:** cad-20562, confirmed housing but retains existing 3.18m bounded plinth.\n\n**Second pilot:** ${secondaryIds.join(', ')}. These seven gray-ridged apartment bars are in a dense western estate, despite the source rural tag; all already seated flat. No blanket rural restoration.\n\n## Image evidence and alignment\n\nSeven local north-up z19 satellite mosaics, each 1792×1792; all 343 tiles decoded, zero missing tiles. Local .png files may contain JPEG; both decoded. Resolution approximately 0.263m/pixel (not an assertion of acquisition resolution). Acquisition date unknown. Green = active surveyed CAD outline; red = audit candidate, not automatic approval. Labels are CAD IDs. All seven outlined images and the two primary unmarked estate images were actually opened using read_image.\n\nThe confirmed slab and industrial roof outlines visually align; roof-to-footprint parallax remains. cad-27730 does **not** match one roof. cad-20346 is a **covered duplicate of another nonactive candidate**, not an active-ID omission to export independently.\n\n${results.map((r) => `- ${r.name}: ${link(`${r.name}-satellite.png`, 'satellite')} · ${link(`${r.name}-outlines.png`, 'CAD overlay')} · ${link(`${r.name}-evidence.json`, 'tile provenance and seating')}`).join('\n')}\n\n## Candidate decisions (20 IDs)\n\n| ID | WGS84 lon, lat | Area m² | Relief m / seating | Decision | Evidence |\n|---|---|---:|---|---|---|\n${rows.join('\n')}\n\n${details.join('\n\n')}\n\n## Near architecture / roof guidance\n\nThe residential pilot roofs are gray, narrow and ridged; maintain per-building direction and actual notched footprint. Use bounded near LOD ridge/eave relief and subtle roof-scale variation, not arbitrary tower floors. Industrial blue roofs have very different sheet/seam morphology and must not receive apartment detailing. Roof-only satellite sampling must mask to its intended surface and avoid baking displaced shadows, trees or roads; never claim satellite-derived facade windows or measured height. Exact floors, facade materials and heights are uncertain.\n\n## Scope and validation\n\nMachine-readable decisions: ${link('review.json')}. Tool: [review-large-building-imagery.mjs](<tools/city-massing/review-large-building-imagery.mjs>). Report-generation assertions enforce already-seated and unique pilot candidates. Other-neighborhood zero overlap was checked; runtime browser appearance was not. No plans, assets/config or converter/assembler files changed. Offline image generation and geometry overlap checks only. **Build/dev/server/browser debug skipped under AGENTS.md.**\n`)
console.log(JSON.stringify({ review: join(out, 'review.json'), reviewed: reviewed.length, primaryIds, extendedIds, secondaryIds }))