City-Lifeline-Cockpit/tools/city-massing/verify-imagery-footprints.mjs

// Full-resolution visual check of imagery-derived footprints: crops z19 imagery 1:1 and outlines
// the accepted footprints so a human can judge whether they really sit on roofs.
import { readFileSync, writeFileSync, existsSync } from 'node:fs'
import { resolve } from 'node:path'
import { PNG } from 'pngjs'
import jpeg from 'jpeg-js'
import * as turf from '@turf/turf'

const DIR = process.env.GIS_DIR || 'D:/GIS/GIS'
const cityDir = resolve('tools/city-massing/city')
const collections = [
    { path: resolve(cityDir, 'imagery-footprints.geojson'), color: [255, 235, 0], label: 'imagery' },
    { path: resolve(cityDir, 'osm-buildings-clipped.geojson'), color: [0, 220, 255], label: 'osm' },
    { path: resolve(cityDir, 'regular-estate-roofs.geojson'), color: [255, 140, 0], label: 'estate' },
    { path: resolve(cityDir, 'roads-clipped.geojson'), color: [255, 0, 255], label: 'roads' },
]
const layers = collections
    .filter((entry) => existsSync(entry.path))
    .map((entry) => ({ ...entry, features: JSON.parse(readFileSync(entry.path)).features }))

const ZOOM = 19
const TILE = 256
const CROP = 4 // tiles per side
const panels = []
const centres = JSON.parse(process.argv[2] || '[]')

function project([lon, lat]) {
    const scale = 2 ** ZOOM
    return [
        ((lon + 180) / 360) * scale * TILE,
        ((1 - Math.asinh(Math.tan((lat * Math.PI) / 180)) / Math.PI) / 2) * scale * TILE,
    ]
}

for (const centre of centres) {
    const [cx, cy] = project(centre)
    const x0 = Math.floor(cx / TILE) - CROP / 2 + 1
    const y0 = Math.floor(cy / TILE) - CROP / 2 + 1
    const size = CROP * TILE
    const image = new PNG({ width: size, height: size })
    image.data.fill(30)
    for (let ty = 0; ty < CROP; ty++) {
        for (let tx = 0; tx < CROP; tx++) {
            const gx = x0 + tx, gy = y0 + ty
            const file = `${DIR}/yuanlingPic3/${ZOOM}/${gx}/${2 ** ZOOM - 1 - gy}.png`
            if (!existsSync(file)) continue
            const bytes = readFileSync(file)
            const decoded = bytes[0] === 0xff ? jpeg.decode(bytes) : PNG.sync.read(bytes)
            PNG.bitblt(decoded, image, 0, 0, TILE, TILE, tx * TILE, ty * TILE)
        }
    }
    const originX = x0 * TILE, originY = y0 * TILE
    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) - originX
            const y = Math.round(a[1] + ((b[1] - a[1]) * i) / steps) - originY
            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 >= size || py >= size) continue
                    const k = (py * size + px) * 4
                    image.data[k] = color[0]
                    image.data[k + 1] = color[1]
                    image.data[k + 2] = color[2]
                    image.data[k + 3] = 255
                }
        }
    }
    const window = turf.bboxPolygon([
        ...turf.bbox(turf.point(centre)).slice(0, 2),
        centre[0] + 0.004,
        centre[1] + 0.004,
    ])
    let drawn = 0
    for (const layer of layers) {
        for (const feature of layer.features) {
            const type = feature.geometry.type
            const rings = type === 'Polygon' ? [feature.geometry.coordinates[0]]
                : type === 'MultiPolygon' ? feature.geometry.coordinates.map((poly) => poly[0])
                    : type === 'LineString' ? [feature.geometry.coordinates]
                        : type === 'MultiLineString' ? feature.geometry.coordinates
                            : []
            if (!rings.length) continue
            if (!turf.booleanIntersects(feature, window)) continue
            for (const ring of rings)
                for (let i = 0; i < ring.length - 1; i++) draw(project(ring[i]), project(ring[i + 1]), layer.color, 1)
            drawn++
        }
    }
    panels.push({ image, drawn, centre })
}

const totalWidth = Math.max(...panels.map((panel) => panel.image.width))
const totalHeight = panels.reduce((sum, panel) => sum + panel.image.height, 0)
const sheet = new PNG({ width: totalWidth, height: totalHeight })
sheet.data.fill(20)
let offset = 0
for (const panel of panels) {
    PNG.bitblt(panel.image, sheet, 0, 0, panel.image.width, panel.image.height, 0, offset)
    offset += panel.image.height
}
const out = resolve(cityDir, 'imagery-footprints-verify.png')
writeFileSync(out, PNG.sync.write(sheet))
console.log(JSON.stringify({
    output: out,
    panels: panels.map((panel) => ({ centre: panel.centre, outlinesDrawn: panel.drawn })),
    legend: 'yellow = imagery footprints, cyan = mapped OSM footprints',
}))