232 lines
11 KiB
JavaScript
232 lines
11 KiB
JavaScript
"use strict";
|
|
|
|
const fs = require("fs");
|
|
|
|
const PI = Math.PI;
|
|
const EARTH_A = 6378245.0;
|
|
const EARTH_EE = 0.00669342162296594323;
|
|
|
|
function transformLat(x, y) {
|
|
let value = -100 + 2 * x + 3 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * Math.sqrt(Math.abs(x));
|
|
value += (20 * Math.sin(6 * x * PI) + 20 * Math.sin(2 * x * PI)) * 2 / 3;
|
|
value += (20 * Math.sin(y * PI) + 40 * Math.sin(y / 3 * PI)) * 2 / 3;
|
|
value += (160 * Math.sin(y / 12 * PI) + 320 * Math.sin(y * PI / 30)) * 2 / 3;
|
|
return value;
|
|
}
|
|
|
|
function transformLon(x, y) {
|
|
let value = 300 + x + 2 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * Math.sqrt(Math.abs(x));
|
|
value += (20 * Math.sin(6 * x * PI) + 20 * Math.sin(2 * x * PI)) * 2 / 3;
|
|
value += (20 * Math.sin(x * PI) + 40 * Math.sin(x / 3 * PI)) * 2 / 3;
|
|
value += (150 * Math.sin(x / 12 * PI) + 300 * Math.sin(x / 30 * PI)) * 2 / 3;
|
|
return value;
|
|
}
|
|
|
|
// This is the standard local inverse approximation used for GCJ-02 reference
|
|
// data. It is intentionally kept separate from native road geometry, whose
|
|
// source coordinates remain WGS84.
|
|
function gcj02ToWgs84(coordinate) {
|
|
const [longitude, latitude] = coordinate;
|
|
if (!Number.isFinite(longitude) || !Number.isFinite(latitude)) throw new Error("Reference coordinate must be finite");
|
|
const dLat = transformLat(longitude - 105, latitude - 35);
|
|
const dLon = transformLon(longitude - 105, latitude - 35);
|
|
const radLat = latitude / 180 * PI;
|
|
const magic = 1 - EARTH_EE * Math.sin(radLat) ** 2;
|
|
const sqrtMagic = Math.sqrt(magic);
|
|
return [
|
|
longitude - dLon * 180 / (EARTH_A / sqrtMagic * Math.cos(radLat) * PI),
|
|
latitude - dLat * 180 / (EARTH_A * (1 - EARTH_EE) / (magic * sqrtMagic) * PI),
|
|
];
|
|
}
|
|
|
|
function mapCoordinates(coordinates, mapper) {
|
|
if (typeof coordinates[0] === "number") return mapper(coordinates);
|
|
return coordinates.map((value) => mapCoordinates(value, mapper));
|
|
}
|
|
|
|
function convertGeoJson(document) {
|
|
if (!document || document.type !== "FeatureCollection" || !Array.isArray(document.features)) {
|
|
throw new Error("Reference must be a GeoJSON FeatureCollection");
|
|
}
|
|
return {
|
|
...document,
|
|
crs: undefined,
|
|
features: document.features.map((feature) => {
|
|
if (!feature || !feature.geometry || !feature.geometry.coordinates) throw new Error("Reference feature is missing geometry");
|
|
return { ...feature, geometry: { ...feature.geometry, coordinates: mapCoordinates(feature.geometry.coordinates, gcj02ToWgs84) } };
|
|
}),
|
|
};
|
|
}
|
|
|
|
function coordinatesOf(document) {
|
|
const points = [];
|
|
for (const feature of document.features || []) walkCoordinates(feature.geometry?.coordinates, points);
|
|
return points;
|
|
}
|
|
|
|
function walkCoordinates(value, points) {
|
|
if (!Array.isArray(value) || !value.length) return;
|
|
if (typeof value[0] === "number") {
|
|
points.push(value);
|
|
return;
|
|
}
|
|
for (const child of value) walkCoordinates(child, points);
|
|
}
|
|
|
|
function boundsOf(document) {
|
|
const points = coordinatesOf(document);
|
|
if (!points.length) throw new Error("Reference contains no coordinates");
|
|
return {
|
|
minLon: Math.min(...points.map((point) => point[0])),
|
|
minLat: Math.min(...points.map((point) => point[1])),
|
|
maxLon: Math.max(...points.map((point) => point[0])),
|
|
maxLat: Math.max(...points.map((point) => point[1])),
|
|
};
|
|
}
|
|
|
|
function centerOf(bounds) {
|
|
return [(bounds.minLon + bounds.maxLon) / 2, (bounds.minLat + bounds.maxLat) / 2];
|
|
}
|
|
|
|
function distanceMeters(first, second) {
|
|
const lonScale = 111320 * Math.cos(first[1] * PI / 180);
|
|
return Math.hypot((second[0] - first[0]) * lonScale, (second[1] - first[1]) * 111320);
|
|
}
|
|
|
|
function parseOsmNodes(xml) {
|
|
const nodes = [];
|
|
for (const match of xml.matchAll(/<node\b([^>]*?)(?:\/>|>([\s\S]*?)<\/node>)/g)) {
|
|
const attrs = {};
|
|
for (const item of match[1].matchAll(/([:\w-]+)\s*=\s*(?:"([^"]*)"|'([^']*)')/g)) attrs[item[1]] = item[2] ?? item[3];
|
|
if (!attrs.id || !Number.isFinite(Number(attrs.lon)) || !Number.isFinite(Number(attrs.lat))) continue;
|
|
const tags = {};
|
|
for (const item of (match[2] || "").matchAll(/<tag\b([^>]*)\/?\s*>/g)) {
|
|
const tag = {};
|
|
for (const attr of item[1].matchAll(/([:\w-]+)\s*=\s*(?:"([^"]*)"|'([^']*)')/g)) tag[attr[1]] = attr[2] ?? attr[3];
|
|
if (tag.k) tags[tag.k] = tag.v || "";
|
|
}
|
|
nodes.push({ id: String(attrs.id), coordinate: [Number(attrs.lon), Number(attrs.lat)], tags });
|
|
}
|
|
return nodes;
|
|
}
|
|
|
|
function nearestNode(nodes, coordinate, nodeId) {
|
|
if (nodeId) {
|
|
const exact = nodes.find((node) => node.id === String(nodeId));
|
|
if (!exact) throw new Error(`OSM node not found: ${nodeId}`);
|
|
return { ...exact, distanceMeters: distanceMeters(exact.coordinate, coordinate), match: "node-id" };
|
|
}
|
|
const candidates = nodes.map((node) => ({ ...node, distanceMeters: distanceMeters(node.coordinate, coordinate) }));
|
|
candidates.sort((first, second) => first.distanceMeters - second.distanceMeters);
|
|
if (!candidates[0]) throw new Error("OSM contains no usable nodes");
|
|
return { ...candidates[0], match: "nearest-node" };
|
|
}
|
|
|
|
function bboxIntersectionRatio(first, second) {
|
|
const width = Math.max(0, Math.min(first.maxLon, second.maxLon) - Math.max(first.minLon, second.minLon));
|
|
const height = Math.max(0, Math.min(first.maxLat, second.maxLat) - Math.max(first.minLat, second.minLat));
|
|
const intersection = width * height;
|
|
const firstArea = Math.max(0, first.maxLon - first.minLon) * Math.max(0, first.maxLat - first.minLat);
|
|
const secondArea = Math.max(0, second.maxLon - second.minLon) * Math.max(0, second.maxLat - second.minLat);
|
|
return intersection / Math.max(firstArea + secondArea - intersection, Number.EPSILON);
|
|
}
|
|
|
|
// A complex junction is compiled as one cluster of `complex_part` polygons in
|
|
// road_surface.geojson, not as a per-node feature in intersection_surface.
|
|
// Match it by cluster id, or by whichever cluster core sits nearest the node.
|
|
function complexClusterSurface(nativeRoadSurfaceFile, node, clusterId) {
|
|
if (!nativeRoadSurfaceFile || !fs.existsSync(nativeRoadSurfaceFile)) return null;
|
|
const surface = JSON.parse(fs.readFileSync(nativeRoadSurfaceFile, "utf8"));
|
|
const parts = (surface.features || []).filter((item) => item.properties?.cluster_id && item.properties?.complex_part);
|
|
const cores = parts.filter((item) => item.properties.complex_part === "core" && Array.isArray(item.properties.center));
|
|
if (!cores.length) return null;
|
|
const core = clusterId
|
|
? cores.find((item) => String(item.properties.cluster_id) === String(clusterId))
|
|
: [...cores].sort((first, second) => distanceMeters(first.properties.center, node.coordinate) - distanceMeters(second.properties.center, node.coordinate))[0];
|
|
if (!core) return null;
|
|
const features = parts.filter((item) => item.properties.cluster_id === core.properties.cluster_id);
|
|
return { clusterId: core.properties.cluster_id, core, features };
|
|
}
|
|
|
|
function inspectReference({ referenceFile, osmFile, nativeIntersectionFile, nativeRoadSurfaceFile, nodeId, clusterId }) {
|
|
const source = JSON.parse(fs.readFileSync(referenceFile, "utf8"));
|
|
const converted = convertGeoJson(source);
|
|
const referenceBounds = boundsOf(converted);
|
|
const referenceCenter = centerOf(referenceBounds);
|
|
const nodes = parseOsmNodes(fs.readFileSync(osmFile, "utf8"));
|
|
const matchedNode = nearestNode(nodes, referenceCenter, nodeId);
|
|
const native = JSON.parse(fs.readFileSync(nativeIntersectionFile, "utf8"));
|
|
const feature = (native.features || []).find((item) => item.properties?.osm_node_id === matchedNode.id);
|
|
const cluster = feature ? null : complexClusterSurface(nativeRoadSurfaceFile, matchedNode, clusterId);
|
|
const matchedFeatures = feature ? [feature] : cluster?.features || null;
|
|
const nativeBounds = matchedFeatures ? boundsOf({ features: matchedFeatures }) : null;
|
|
const diagnostics = [];
|
|
if (!matchedFeatures) diagnostics.push(nativeRoadSurfaceFile ? "No native intersection surface or complex cluster matched the OSM node" : "No native intersection surface matched the OSM node; pass --native-road-surface to also search complex junction clusters");
|
|
return {
|
|
schema: "gaode-junction-reference-comparison/v2",
|
|
source: { file: referenceFile, coordinateSystem: "GCJ-02", featureCount: converted.features.length },
|
|
conversion: { target: "WGS84", method: "gcj02-inverse-approximation" },
|
|
reference: { bounds: referenceBounds, center: referenceCenter },
|
|
matchedOsmNode: { id: matchedNode.id, coordinate: matchedNode.coordinate, tags: matchedNode.tags, match: matchedNode.match, centerDistanceMeters: matchedNode.distanceMeters },
|
|
nativeIntersection: nativeBounds ? {
|
|
kind: feature ? "junction-node" : "complex-cluster",
|
|
clusterId: cluster?.clusterId || null,
|
|
featureCount: matchedFeatures.length,
|
|
bounds: nativeBounds,
|
|
bboxIoU: bboxIntersectionRatio(referenceBounds, nativeBounds),
|
|
centerOffsetMeters: distanceMeters(referenceCenter, centerOf(nativeBounds)),
|
|
featureProperties: feature ? feature.properties : cluster.core.properties,
|
|
} : null,
|
|
diagnostics,
|
|
converted,
|
|
matchedFeatures,
|
|
};
|
|
}
|
|
|
|
function localReferenceSvg({ converted, nativeIntersection, center, radiusMeters = 180 }) {
|
|
const width = 1000;
|
|
const height = 1000;
|
|
const lonScale = 111320 * Math.cos(center[1] * PI / 180);
|
|
const project = (point) => [
|
|
width / 2 + (point[0] - center[0]) * lonScale * width / (radiusMeters * 2),
|
|
height / 2 - (point[1] - center[1]) * 111320 * height / (radiusMeters * 2),
|
|
];
|
|
const pathFor = (coordinates) => {
|
|
const parts = [];
|
|
const appendLine = (line, close) => {
|
|
if (!line?.length) return;
|
|
const [firstX, firstY] = project(line[0]);
|
|
parts.push(`M ${firstX.toFixed(1)} ${firstY.toFixed(1)}`);
|
|
for (const point of line.slice(1)) {
|
|
const [x, y] = project(point);
|
|
parts.push(`L ${x.toFixed(1)} ${y.toFixed(1)}`);
|
|
}
|
|
if (close) parts.push("Z");
|
|
};
|
|
const visit = (value) => {
|
|
if (!Array.isArray(value) || !value.length) return;
|
|
if (typeof value[0] === "number") return;
|
|
if (typeof value[0][0] === "number") appendLine(value, value.length > 2);
|
|
else value.forEach(visit);
|
|
};
|
|
visit(coordinates);
|
|
return parts.join(" ");
|
|
};
|
|
const color = { 1: "#2563eb", 2: "#0f766e", 3: "#7c3aed", 4: "#ea580c", 5: "#64748b" };
|
|
const references = converted.features.map((feature) => {
|
|
const type = feature.properties?.type || "unknown";
|
|
return `<path d="${pathFor(feature.geometry.coordinates)}" fill="${feature.geometry.type.includes("Polygon") ? `${color[type] || "#334155"}18` : "none"}" stroke="${color[type] || "#334155"}" stroke-width="1.2"/>`;
|
|
}).join("\n");
|
|
const nativePaths = (nativeIntersection?.features || []).map((feature) => `<path d="${pathFor(feature.geometry.coordinates)}" fill="#dc262655" stroke="#dc2626" stroke-width="3"/>`).join("\n");
|
|
return `<?xml version="1.0" encoding="UTF-8"?>
|
|
<svg xmlns="http://www.w3.org/2000/svg" width="${width}" height="${height}" viewBox="0 0 ${width} ${height}">
|
|
<rect width="100%" height="100%" fill="#f8fafc"/>
|
|
${references}
|
|
${nativePaths}
|
|
<circle cx="${width / 2}" cy="${height / 2}" r="5" fill="#111827"/>
|
|
<text x="20" y="35" font-family="sans-serif" font-size="22" fill="#111827">Gaode reference (type colors) / native intersection (red)</text>
|
|
</svg>`;
|
|
}
|
|
|
|
module.exports = { gcj02ToWgs84, convertGeoJson, boundsOf, parseOsmNodes, nearestNode, inspectReference, localReferenceSvg };
|