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