// 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 }))