City-Lifeline-Cockpit/tools/city-massing/city-atlas.mjs

// Offline city inventory. WGS84 bounds are the explicitly agreed deployment area, not GIS tile coverage.
import { readFileSync, writeFileSync, mkdirSync, existsSync } from 'node:fs'
import { resolve } from 'node:path'
import { PNG } from 'pngjs'
import jpeg from 'jpeg-js'
import * as turf from '@turf/turf'
import { PRIMARY_AREA_BOUNDS, clipFeatureToBounds } from './primary-area.mjs'

const DIR = process.env.GIS_DIR || 'D:/GIS/GIS'
const OUT = resolve('tools/city-massing/city')
const BOUNDS = PRIMARY_AREA_BOUNDS
const ZOOM = 15
const TILE = 256
mkdirSync(OUT, { recursive: true })
const slippy = ([lon, lat]) => [
    ((lon + 180) / 360) * 2 ** ZOOM,
    ((1 - Math.asinh(Math.tan((lat * Math.PI) / 180)) / Math.PI) / 2) * 2 ** ZOOM,
]
const nw = slippy([BOUNDS[0], BOUNDS[3]])
const se = slippy([BOUNDS[2], BOUNDS[1]])
const [west, north, east, south] = [Math.floor(nw[0]), Math.floor(nw[1]), Math.ceil(se[0]), Math.ceil(se[1])]
const image = new PNG({ width: (east - west) * TILE, height: (south - north) * TILE })
image.data.fill(230)
let tiles = 0
const missing = []
for (let y = north; y < south; y++)
    for (let x = west; x < east; x++) {
        const file = `${DIR}/yuanlingPic3/${ZOOM}/${x}/${2 ** ZOOM - 1 - y}.png`
        if (!existsSync(file)) {
            missing.push([x, y])
            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, (x - west) * TILE, (y - north) * TILE)
        tiles++
    }
if (missing.length) throw new Error(`City atlas incomplete: ${missing.length} image tiles missing`)
writeFileSync(resolve(OUT, 'city-satellite.png'), PNG.sync.write(image))
const window = turf.bboxPolygon(BOUNDS)
function clippedFeatures(source) {
    const featureCollection = JSON.parse(readFileSync(source, 'utf8'))
    const result = []
    for (const feature of featureCollection.features) {
        try {
            const clipped = clipFeatureToBounds(feature, BOUNDS)
            if (clipped) result.push(clipped)
        } catch (error) {
            throw new Error(`Cannot clip ${source}: ${error.message}`)
        }
    }
    return result
}
const residential = clippedFeatures(`${DIR}/output/resa.json`)
const water = clippedFeatures(`${DIR}/output/hyda.json`)
const osmBase = `${DIR}/planet_110.059,28.258_110.744,28.625-shp/shpOutput`
const osmBuildings = clippedFeatures(`${osmBase}/buildings.json`)
const osmRoads = clippedFeatures(`${osmBase}/roads.json`)
const osmLanduse = clippedFeatures(`${osmBase}/landuse.json`)
const sections = [
    { id: 'south', label: '城南', bbox: [BOUNDS[0], BOUNDS[1], BOUNDS[2], 28.432] },
    { id: 'river-central', label: '沿江老城区', bbox: [BOUNDS[0], 28.432, BOUNDS[2], 28.455] },
    { id: 'north', label: '城北', bbox: [BOUNDS[0], 28.455, BOUNDS[2], BOUNDS[3]] },
]
const districts = sections.map((section) => {
    const bbox = turf.bboxPolygon(section.bbox)
    const inside = (xs) => xs.filter((x) => turf.booleanIntersects(x, bbox))
    const masks = inside(residential)
    return {
        ...section,
        residentialIds: masks.map((x) => x.properties.ORIG_FID),
        residentialSqKm: Math.round(masks.reduce((s, x) => s + turf.area(x), 0) / 10000) / 100,
        osmBuildingCount: inside(osmBuildings).length,
        roadCount: inside(osmRoads).length,
        status: 'unreviewed',
        priority: section.id === 'south' ? 'regular-first' : 'old-town-and-regular',
    }
})
const atlas = {
    source: 'yuanlingPic3 TMS + output/resa.json residential polygon + OSM buildings/roads',
    bounds: BOUNDS,
    imagery: {
        zoom: ZOOM,
        westTile: west,
        googleY0: north,
        tmsTopY: 2 ** ZOOM - 1 - north,
        tileCount: tiles,
        width: image.width,
        height: image.height,
    },
    coveragePolicy:
        'Building coverage is measured only within independently verified urban masks; this inventory is not itself a modeled-building footprint catalog.',
    districts,
    originalSources: {
        residential: residential.length,
        osmBuildings: osmBuildings.length,
        osmRoads: osmRoads.length,
        water: water.length,
        osmLanduse: osmLanduse.length,
    },
}
writeFileSync(resolve(OUT, 'city-inventory.json'), JSON.stringify(atlas, null, 2))
writeFileSync(resolve(OUT, 'residential-clipped.geojson'), JSON.stringify(turf.featureCollection(residential)))
writeFileSync(resolve(OUT, 'water-clipped.geojson'), JSON.stringify(turf.featureCollection(water)))
writeFileSync(resolve(OUT, 'osm-buildings-clipped.geojson'), JSON.stringify(turf.featureCollection(osmBuildings)))
writeFileSync(resolve(OUT, 'roads-clipped.geojson'), JSON.stringify(turf.featureCollection(osmRoads)))
console.log(JSON.stringify({ tiles, atlasPixels: [image.width, image.height], ...atlas.originalSources, districts }))