City-Lifeline-Cockpit/tools/city-massing/draw-city-masks.mjs

// Draw source residential masks, waterways and main roads over the georeferenced city atlas.
import { readFileSync, writeFileSync } from 'node:fs'
import { resolve } from 'node:path'
import { PNG } from 'pngjs'
const dir = resolve('tools/city-massing/city')
const atlas = JSON.parse(readFileSync(resolve(dir, 'city-inventory.json')))
const image = PNG.sync.read(readFileSync(resolve(dir, 'city-satellite.png')))
const b = atlas.bounds,
    z = atlas.imagery.zoom,
    n = 2 ** z,
    tile = 256
function project([lon, lat]) {
    const x = ((lon + 180) / 360) * n * tile,
        y = ((1 - Math.asinh(Math.tan((lat * Math.PI) / 180)) / Math.PI) / 2) * n * tile
    return [x - atlas.imagery.westTile * tile, y - atlas.imagery.googleY0 * tile]
}
function line(a, b, color, width = 2) {
    const len = Math.ceil(Math.hypot(b[0] - a[0], b[1] - a[1]))
    for (let i = 0; i <= len; i++) {
        const x = Math.round(a[0] + ((b[0] - a[0]) * i) / len),
            y = Math.round(a[1] + ((b[1] - a[1]) * i) / len)
        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
            }
    }
}
const layers = [
    ['residential-clipped.geojson', [255, 220, 0], 3],
    ['water-clipped.geojson', [0, 210, 255], 3],
    ['osm-buildings-clipped.geojson', [255, 30, 30], 2],
]
for (const [file, color, width] of layers) {
    const j = JSON.parse(readFileSync(resolve(dir, file)))
    for (const f of j.features) {
        if (f.geometry.type === 'Polygon')
            for (const ring of f.geometry.coordinates)
                for (let i = 0; i < ring.length - 1; i++) line(project(ring[i]), project(ring[i + 1]), color, width)
        if (f.geometry.type === 'MultiPolygon')
            for (const poly of f.geometry.coordinates)
                for (const ring of poly)
                    for (let i = 0; i < ring.length - 1; i++) line(project(ring[i]), project(ring[i + 1]), color, width)
    }
}
const roads = JSON.parse(readFileSync(resolve(dir, 'roads-clipped.geojson')))
for (const f of roads.features) {
    const c = f.geometry.coordinates
    if (f.geometry.type === 'LineString')
        for (let i = 1; i < c.length; i++) line(project(c[i - 1]), project(c[i]), [255, 0, 255], 1)
    else if (f.geometry.type === 'MultiLineString')
        for (const l of c) for (let i = 1; i < l.length; i++) line(project(l[i - 1]), project(l[i]), [255, 0, 255], 1)
}
writeFileSync(resolve(dir, 'city-inventory-overlay.png'), PNG.sync.write(image))
console.log(
    JSON.stringify({
        size: [image.width, image.height],
        bbox: b,
        imagery: atlas.imagery,
        layerCounts: {
            residential: atlas.originalSources.residential,
            water: atlas.originalSources.water,
            roads: atlas.originalSources.osmRoads,
            osmBuildings: atlas.originalSources.osmBuildings,
        },
    })
)