City-Lifeline-Cockpit/tools/city-massing/verify-cad-buildings.mjs

// Validate the CAD building layer against the local imagery and report what the survey contains.
// Writes an overlay so the alignment claim can be checked by eye, not just by counts.
import { readFileSync, writeFileSync } from 'node:fs'
import { resolve } from 'node:path'
import { PNG } from 'pngjs'
import * as turf from '@turf/turf'
import { PRIMARY_AREA_BOUNDS, isWithinPrimaryArea } from './primary-area.mjs'

const cityDir = resolve('tools/city-massing/city')
const buildings = JSON.parse(readFileSync(resolve(cityDir, 'cad-buildings.geojson'), 'utf8')).features
const labelsPath = resolve(cityDir, 'cad-building-labels.geojson')
let labels = []
try {
    labels = JSON.parse(readFileSync(labelsPath, 'utf8')).features
} catch {
    labels = []
}

const areas = []
let outside = 0
let degenerate = 0
let vertexHistogram = {}
for (const feature of buildings) {
    const ring = feature.geometry.coordinates[0]
    vertexHistogram[ring.length] = (vertexHistogram[ring.length] || 0) + 1
    const area = turf.area(feature)
    if (!(area > 0)) {
        degenerate++
        continue
    }
    areas.push(area)
    if (!isWithinPrimaryArea(feature)) outside++
}
areas.sort((a, b) => a - b)
const pick = (q) => Math.round(areas[Math.floor(q * (areas.length - 1))] * 10) / 10

// Label text: survey annotations often carry storeys and structure (e.g. "砖混3").
const counts = {}
let withStructure = 0
let withStoreys = 0
for (const label of labels) {
    const text = (label.properties.text || '').trim()
    if (!text) continue
    counts[text] = (counts[text] || 0) + 1
    if (/砖混|砖木|框架|土木|钢/.test(text)) withStructure++
    if (/\d/.test(text)) withStoreys++
}
const topLabels = Object.entries(counts).sort((a, b) => b[1] - a[1]).slice(0, 25)

// --- overlay on the city atlas -------------------------------------------------------------
const atlas = JSON.parse(readFileSync(resolve(cityDir, 'city-inventory.json')))
const image = PNG.sync.read(readFileSync(resolve(cityDir, 'city-satellite.png')))
const scale = 2 ** atlas.imagery.zoom * 256
const project = ([lon, lat]) => [
    ((lon + 180) / 360) * scale - atlas.imagery.westTile * 256,
    ((1 - Math.asinh(Math.tan((lat * Math.PI) / 180)) / Math.PI) / 2) * scale - atlas.imagery.googleY0 * 256,
]
function draw(a, b, color, width) {
    const steps = Math.max(1, 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)
        const y = Math.round(a[1] + ((b[1] - a[1]) * i) / steps)
        for (let dy = -width; dy <= width; dy++)
            for (let dx = -width; dx <= width; dx++) {
                const px = x + dx, py = y + dy
                if (px < 0 || py < 0 || px >= image.width || py >= image.height) continue
                const k = (py * image.width + px) * 4
                image.data[k] = color[0]
                image.data[k + 1] = color[1]
                image.data[k + 2] = color[2]
                image.data[k + 3] = 255
            }
    }
}
let drawn = 0
for (const feature of buildings) {
    const rings = feature.geometry.type === 'Polygon' ? feature.geometry.coordinates : []
    for (const ring of rings) {
        const points = ring.map(project)
        for (let i = 0; i < points.length - 1; i++) draw(points[i], points[i + 1], [255, 60, 200], 0)
    }
    drawn++
}
writeFileSync(resolve(cityDir, 'cad-buildings-overlay.png'), PNG.sync.write(image))
// 1:1 z19 crops are the only way to judge whether a footprint really sits on its roof.
import { readFileSync as read19, existsSync as exists19 } from 'node:fs'
import jpeg19 from 'jpeg-js'
const Z19 = 19, TILE19 = 256, CROP_TILES = 4
function project19([lon, lat]) {
    const s = 2 ** Z19
    return [((lon + 180) / 360) * s * TILE19, ((1 - Math.asinh(Math.tan((lat * Math.PI) / 180)) / Math.PI) / 2) * s * TILE19]
}
function crop19(centre, outName) {
    const [cx19, cy19] = project19(centre)
    const x0 = Math.floor(cx19 / TILE19) - CROP_TILES / 2 + 1
    const y0 = Math.floor(cy19 / TILE19) - CROP_TILES / 2 + 1
    const size = CROP_TILES * TILE19
    const canvas = new PNG({ width: size, height: size })
    canvas.data.fill(30)
    const dir = process.env.GIS_DIR || 'D:/GIS/GIS'
    for (let ty = 0; ty < CROP_TILES; ty++) {
        for (let tx = 0; tx < CROP_TILES; tx++) {
            const gx = x0 + tx, gy = y0 + ty
            const file = `${dir}/yuanlingPic3/${Z19}/${gx}/${2 ** Z19 - 1 - gy}.png`
            if (!exists19(file)) continue
            const bytes = read19(file)
            const decoded = bytes[0] === 0xff ? jpeg19.decode(bytes) : PNG.sync.read(bytes)
            PNG.bitblt(decoded, canvas, 0, 0, TILE19, TILE19, tx * TILE19, ty * TILE19)
        }
    }
    const originX = x0 * TILE19, originY = y0 * TILE19
    let n = 0
    for (const feature of buildings) {
        for (const ring of feature.geometry.type === 'Polygon' ? feature.geometry.coordinates : []) {
            const pts = ring.map((p) => project19(p))
            if (!pts.some(([px, py]) => px > originX - 256 && px < originX + size + 256 && py > originY - 256 && py < originY + size + 256)) continue
            n++
            for (let i = 0; i < pts.length - 1; i++) {
                const a = [pts[i][0] - originX, pts[i][1] - originY]
                const b = [pts[i + 1][0] - originX, pts[i + 1][1] - originY]
                const steps = Math.max(1, Math.ceil(Math.hypot(b[0] - a[0], b[1] - a[1])))
                for (let k = 0; k <= steps; k++) {
                    const x = Math.round(a[0] + ((b[0] - a[0]) * k) / steps)
                    const y = Math.round(a[1] + ((b[1] - a[1]) * k) / steps)
                    if (x < 0 || y < 0 || x >= size || y >= size) continue
                    const off = (y * size + x) * 4
                    canvas.data[off] = 255; canvas.data[off + 1] = 40; canvas.data[off + 2] = 200; canvas.data[off + 3] = 255
                }
            }
        }
    }
    writeFileSync(resolve(cityDir, outName), PNG.sync.write(canvas))
    return { out: outName, outlines: n }
}
const crops = [crop19([110.3965, 28.4575], 'cad-buildings-z19-core.png')]

const CROP = 1400
const [cx, cy] = project([110.3965, 28.4570])
const crop = new PNG({ width: CROP, height: CROP })
PNG.bitblt(image, crop, Math.round(cx - CROP / 2), Math.round(cy - CROP / 2), CROP, CROP, 0, 0)
writeFileSync(resolve(cityDir, 'cad-buildings-core.png'), PNG.sync.write(crop))

console.log(JSON.stringify({
    cadBuildings: buildings.length,
    withValidArea: areas.length,
    degenerate,
    outsidePrimaryArea: outside,
    areaSqm: { min: pick(0), p25: pick(0.25), median: pick(0.5), p75: pick(0.75), p95: pick(0.95), max: pick(1) },
    under15Sqm: areas.filter((a) => a < 15).length,
    over2000Sqm: areas.filter((a) => a > 2000).length,
    vertexHistogram,
    labels: { total: labels.length, withStructureText: withStructure, withDigits: withStoreys, top: topLabels },
    oneToOneCrops: crops,
    overlays: ['cad-buildings-overlay.png', 'cad-buildings-core.png'],
    primaryAreaBounds: [...PRIMARY_AREA_BOUNDS],
}, null, 1))