City-Lifeline-Cockpit/tools/city-massing/trace-regular-estate.mjs

// Georeference individually traced roof axes from the verified TMS mosaic and retain the satellite comparison.
import { readFileSync, writeFileSync } from 'node:fs'
import { resolve } from 'node:path'
import jpeg from 'jpeg-js'
import { PNG } from 'pngjs'
import * as turf from '@turf/turf'

const dir = resolve('tools/city-massing/city')
const config = JSON.parse(readFileSync(resolve(dir, 'regular-estate-roofs.json'), 'utf8'))
const [zoom, x0, y0, size] = [config.zoom, config.tileX0, config.googleY0, config.tileSize]
const pixelToGeo = ([x, y]) => {
    const scale = 2 ** zoom
    const lon = (360 * (x0 + x / size)) / scale - 180
    const lat = (Math.atan(Math.sinh(Math.PI * (1 - (2 * (y0 + y / size)) / scale))) * 180) / Math.PI
    return [lon, lat]
}
const imagePath = resolve('.scratch/gis-investigation-subagent/regular-z19.jpg')
const source = jpeg.decode(readFileSync(imagePath))
const annotated = new PNG({ width: source.width, height: source.height })
annotated.data.set(source.data)
const features = []
function line(a, b, color) {
    const count = Math.ceil(Math.hypot(a[0] - b[0], a[1] - b[1]))
    for (let i = 0; i <= count; i++) {
        const x = Math.round(a[0] + ((b[0] - a[0]) * i) / count)
        const y = Math.round(a[1] + ((b[1] - a[1]) * i) / count)
        if (x < 0 || y < 0 || x >= annotated.width || y >= annotated.height) continue
        const offset = (y * annotated.width + x) * 4
        for (let c = 0; c < 3; c++) annotated.data[offset + c] = color[c]
        annotated.data[offset + 3] = 255
    }
}
for (const roof of config.roofStrips) {
    const [a, b] = roof.axisPx
    const d = Math.hypot(b[0] - a[0], b[1] - a[1])
    if (d < roof.roofWidthPx * 2) throw new Error(`Roof strip too short: ${roof.id}`)
    const perpendicular = [-(b[1] - a[1]) / d, (b[0] - a[0]) / d]
    const [dx, dy] = perpendicular.map((v) => (v * roof.roofWidthPx) / 2)
    const cornersPx = [
        [a[0] + dx, a[1] + dy],
        [b[0] + dx, b[1] + dy],
        [b[0] - dx, b[1] - dy],
        [a[0] - dx, a[1] - dy],
    ]
    const geo = cornersPx.map(pixelToGeo)
    const ring = [...geo, geo[0]]
    for (let j = 0; j < 4; j++) line(cornersPx[j], cornersPx[(j + 1) % 4], [255, 208, 0])
    features.push(
        turf.polygon([ring], {
            id: `regular-${roof.id}`,
            source: 'traced-roof',
            atlas: 'regular-estate-roofs.json',
            floors: roof.floors,
            heightSource: 'estimated',
            archetype: 'slab',
            cornersPx,
            roofAreaSqm: Math.round(turf.area(turf.polygon([ring]))),
        })
    )
}
const overlaps = []
for (let i = 0; i < features.length; i++)
    for (let j = i + 1; j < features.length; j++) {
        if (
            turf.booleanOverlap(features[i], features[j]) ||
            turf.booleanContains(features[i], features[j]) ||
            turf.booleanContains(features[j], features[i])
        )
            overlaps.push([features[i].properties.id, features[j].properties.id])
    }
writeFileSync(resolve(dir, 'regular-estate-audit.png'), PNG.sync.write(annotated))
if (overlaps.length) throw new Error(`Roof traces overlap: ${JSON.stringify(overlaps)}`)
writeFileSync(resolve(dir, 'regular-estate-roofs.geojson'), JSON.stringify(turf.featureCollection(features), null, 2))
console.log(
    JSON.stringify({
        site: config.site,
        roofs: features.length,
        areaSqm: features.reduce((sum, f) => sum + f.properties.roofAreaSqm, 0),
        overlaps,
        photo: imagePath,
    })
)