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