/** * 沅陵轻量体块的独立回验(票 34)。 * * 这套几何是**程序化生成**的,所以"看起来对"不够 —— 这里验的是三件可证伪的事: * 1. **地理配准**:取一个真实 OSM 轮廓的质心,按生成时的同一条投影链推出它在 GLB 里的应有位置, * 再在 GLB 顶点里找得到它(容差 1m)→ 证明"经纬度 → 局部 ENU → glTF Y-up"这条链没写反; * 2. **tileset transform 自洽**:把 GLB 包围盒中心经 transform 送到经纬度,应落在目标城区窗口内; * 3. **GLB 结构合法**:header/JSON/BIN 长度自洽、accessor 与 bufferView 对得上、无 NaN。 * * 用法:node tools/city-massing/verify-town.mjs */ import { readFileSync } from 'node:fs' import { join } from 'node:path' const DATA_DIR = process.env.CITY_DATA_DIR || 'D:/GIS/GIS/cityMassing' const ORIGIN = { lon: 110.39028164, lat: 28.4436787 } const DEG = Math.PI / 180 const M_PER_DEG_LAT = 110540 const M_PER_DEG_LON = 111320 * Math.cos(ORIGIN.lat * DEG) const OSM = 'D:/GIS/GIS/planet_110.059,28.258_110.744,28.625-shp/shpOutput' const fails = [] function check(name, pass, evidence) { console.log(`${pass ? 'PASS' : 'FAIL'} ${name}${evidence ? ` — ${evidence}` : ''}`) if (!pass) fails.push(name) } // ---------------------------------------------------------------- 解析 GLB const buf = readFileSync(join(DATA_DIR, 'town.glb')) const magic = buf.toString('ascii', 0, 4) const version = buf.readUInt32LE(4) const totalLen = buf.readUInt32LE(8) check('GLB-header', magic === 'glTF' && version === 2 && totalLen === buf.length, `${magic} v${version} len=${totalLen}/${buf.length}`) const jsonLen = buf.readUInt32LE(12) const jsonType = buf.toString('ascii', 16, 20) const gltf = JSON.parse(buf.toString('utf8', 20, 20 + jsonLen)) check('GLB-json-chunk', jsonType === 'JSON' && !!gltf.meshes, `${jsonType}, primitives=${gltf.meshes[0].primitives.length}`) const binStart = 20 + jsonLen + 8 const binChunkLen = buf.readUInt32LE(20 + jsonLen) check('GLB-bin-chunk', binStart + binChunkLen <= buf.length, `bin=${binChunkLen}`) function readAccessor(index) { const acc = gltf.accessors[index] const view = gltf.bufferViews[acc.bufferView] const offset = binStart + (view.byteOffset || 0) const out = [] if (acc.componentType === 5126) { for (let i = 0; i < acc.count * 3; i++) out.push(buf.readFloatLE(offset + i * 4)) } else if (acc.componentType === 5125) { for (let i = 0; i < acc.count; i++) out.push(buf.readUInt32LE(offset + i * 4)) } return out } // ---------------------------------------------------------------- 3. 结构 + NaN let vertexCount = 0, triCount = 0, nan = 0 for (const prim of gltf.meshes[0].primitives) { const pos = readAccessor(prim.attributes.POSITION) const idx = readAccessor(prim.indices) vertexCount += pos.length / 3 triCount += idx.length / 3 for (const v of pos) if (!Number.isFinite(v)) nan += 1 const acc = gltf.accessors[prim.attributes.POSITION] check( `POSITION-accessor-minmax-${prim.material}`, Array.isArray(acc.min) && Array.isArray(acc.max) && acc.min.every((v, i) => v <= acc.max[i]), `min=[${acc.min.map((v) => v.toFixed(0))}] max=[${acc.max.map((v) => v.toFixed(0))}]`, ) } check('geometry-sane', vertexCount > 1000 && triCount > 1000 && nan === 0, `顶点 ${vertexCount},三角 ${triCount},NaN ${nan}`) // ---------------------------------------------------------------- 1. 真实轮廓回验 // 生成器把体块**居中**放在轮廓外接盒中心上,所以顶点落在**外接盒的四角**(真实轮廓不做旋转)。 // 因此判据是「GLB 里存在一个顶点落在该轮廓的外接盒角点上」——这条能证伪整条坐标链 // (经纬度 → 局部 ENU → glTF Y-up),写反任何一步都对不上。 // 同时按生成器的规则剔掉贴水的轮廓(那些本来就被避让掉了)。 const buildings = JSON.parse(readFileSync(join(OSM, 'buildings.json'), 'utf8')).features const waterways = JSON.parse(readFileSync(join(OSM, 'waterways.json'), 'utf8')).features const inWindow = (lon, lat) => Math.abs(lon - ORIGIN.lon) <= 0.02 && Math.abs(lat - ORIGIN.lat) <= 0.015 const toLocal = (lon, lat) => ({ x: (lon - ORIGIN.lon) * M_PER_DEG_LON, y: (lat - ORIGIN.lat) * M_PER_DEG_LAT }) const waterLines = [] for (const f of waterways) { const ring = f.geometry && f.geometry.type === 'LineString' ? f.geometry.coordinates : null if (!ring) continue const seg = [] for (const [lon, lat] of ring) { if (inWindow(lon, lat)) seg.push(toLocal(lon, lat)) else if (seg.length) { waterLines.push(seg); seg.length = 0 } } if (seg.length) waterLines.push(seg) } function distanceToWater(p) { let best = Infinity for (const line of waterLines) { for (let i = 1; i < line.length; i++) { const a = line[i - 1], b = line[i] const dx = b.x - a.x, dy = b.y - a.y const len2 = dx * dx + dy * dy || 1 let t = ((p.x - a.x) * dx + (p.y - a.y) * dy) / len2 t = Math.max(0, Math.min(1, t)) best = Math.min(best, Math.hypot(p.x - (a.x + t * dx), p.y - (a.y + t * dy))) } } return best } const targets = [] for (const f of buildings) { const ring = f.geometry && f.geometry.type === 'Polygon' ? f.geometry.coordinates[0] : null if (!ring) continue const pts = ring.filter(([lon, lat]) => inWindow(lon, lat)) if (pts.length < 4) continue const ex = pts.map(([lon]) => (lon - ORIGIN.lon) * M_PER_DEG_LON) const ny = pts.map(([, lat]) => (lat - ORIGIN.lat) * M_PER_DEG_LAT) const minX = Math.min(...ex), maxX = Math.max(...ex) const minY = Math.min(...ny), maxY = Math.max(...ny) const cx = (minX + maxX) / 2, cy = (minY + maxY) / 2 if (distanceToWater({ x: cx, y: cy }) < 22) continue // 生成器会避让贴水的 targets.push({ corners: [[minX, minY], [maxX, minY], [minX, maxY], [maxX, maxY]], area: (maxX - minX) * (maxY - minY) }) } targets.sort((a, b) => b.area - a.area) const sample = targets.slice(0, 8) // 把所有顶点按 ENU 还原(glTF Y-up → ENU:east=x, up=y, north=-z) const allVerts = [] for (const prim of gltf.meshes[0].primitives) { const pos = readAccessor(prim.attributes.POSITION) for (let i = 0; i < pos.length; i += 3) { allVerts.push({ x: pos[i], z: -pos[i + 2] }) // east, north(忽略 up) } } let hit = 0 const details = [] for (const t of sample) { const found = t.corners.some(([cx2, cy2]) => allVerts.some((v) => Math.abs(v.x - cx2) <= 2 && Math.abs(v.z - cy2) <= 2)) if (found) hit += 1 details.push(`${t.area.toFixed(0)}m²@(${t.corners[0][0].toFixed(0)},${t.corners[0][1].toFixed(0)})${found ? '✓' : '✗'}`) } check( 'real-footprints-land-on-true-coords', sample.length >= 5 && hit === sample.length, `${hit}/${sample.length} 个最大真实轮廓的外接盒角点在 GLB 里 2m 内命中:${details.join(' ')}`, ) // ---------------------------------------------------------------- 2. transform 自洽 const T = JSON.parse(readFileSync(join(DATA_DIR, 'tileset.json'), 'utf8')).root.transform const box = JSON.parse(readFileSync(join(DATA_DIR, 'tileset.json'), 'utf8')).root.boundingVolume.box function applyTransform([x, y, z]) { return [ T[0] * x + T[4] * y + T[8] * z + T[12], T[1] * x + T[5] * y + T[9] * z + T[13], T[2] * x + T[6] * y + T[10] * z + T[14], ] } function ecefToLonLat([X, Y, Z]) { const a = 6378137, f = 1 / 298.257223563, e2 = 2 * f - f * f const lon = Math.atan2(Y, X) const p = Math.hypot(X, Y) let lat = Math.atan2(Z, p * (1 - e2)) for (let i = 0; i < 5; i++) { const s = Math.sin(lat) const N = a / Math.sqrt(1 - e2 * s * s) lat = Math.atan2(Z + e2 * N * s, p) } return { lon: (lon * 180) / Math.PI, lat: (lat * 180) / Math.PI } } // 原点自身(局部 0,0,0)应回到目标坐标 const originLL = ecefToLonLat(applyTransform([0, 0, 0])) const dOrigin = Math.hypot((originLL.lon - ORIGIN.lon) * M_PER_DEG_LON, (originLL.lat - ORIGIN.lat) * M_PER_DEG_LAT) check('transform-origin-roundtrip', dOrigin < 0.5, `原点误差 ${dOrigin.toFixed(3)} m(${originLL.lon.toFixed(6)}, ${originLL.lat.toFixed(6)})`) // 包围盒中心经 transform 后应落在窗口内 const corner = ecefToLonLat(applyTransform([box[0], box[1], box[2]])) check( 'bbox-center-in-window', Math.abs(corner.lon - ORIGIN.lon) <= 0.025 && Math.abs(corner.lat - ORIGIN.lat) <= 0.02, `包围盒中心 → ${corner.lon.toFixed(5)}, ${corner.lat.toFixed(5)}`, ) const meta = JSON.parse(readFileSync(join(DATA_DIR, 'meta.json'), 'utf8')) check('disclaimer-present', typeof meta._disclaimer === 'string' && meta._disclaimer.includes('不是实测建筑'), meta._disclaimer.slice(0, 40) + '…') console.log('') if (fails.length) { console.log(`回验未通过:${fails.join(', ')}`) process.exit(1) } console.log(`回验全部通过:${meta.counts.total} 栋体块(真实轮廓 ${meta.counts.real})、${(meta.glbBytes / 1024).toFixed(0)} KB、${meta.materials} 种材质。`)