Resources Numerical accuracy Calculation policy
Builds and reads the map conversion that places a model on a map grid, and re-expresses it in another coordinate system.
Map coordinates go through proj4 in 64-bit. Moving a georeference to another system finds the new rotation and scale by transforming a second point 10 m along the model's X axis, so those two come from a finite difference. Stored eastings, northings and heights are rounded to six decimals of the map unit and axis components to twelve.
6 calculations in 1 source file. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.
Tests tools/bimcraft-test/georef-test.mjs
Rounds an easting, northing or height before it is stored.
round(v × 10⁶) ÷ 10⁶
Six decimals of the map unit: 0.001 mm when the unit is the meter. The next line, r12, does the same to twelve decimals for axis components.
const r6 = (v) => Math.round(v * 1e6) / 1e6;The longitude and latitude of an easting and northing in the chosen system.
The pair is handed to the transformer crs.js builds on proj4. A system already in degrees is returned unchanged.
// ------------------------------------------------------------------ a model with no map conversion (item 1)
export const lonLatOf = (crs, e, n) => (crs.unit === 'deg' ? [e, n] : transformer(crs, null)(e, n));The matrix that takes model coordinates to map coordinates.
x' = s (a x − b y) + E; y' = s (b x + a y) + N; z' = s z + H, with (a, b) the normalized X axis
The axis is normalized with Math.hypot (length 1 is used when it is zero) and multiplied by the scale in 64-bit. Eastings, northings and height go in unchanged.
/** model units → map units (cc-ifc's georef matrix, column-major). */
export function conversionMatrix(conv) {
let ax = conv.XAxisAbscissa ?? 1, ay = conv.XAxisOrdinate ?? 0;
const l = Math.hypot(ax, ay) || 1; ax /= l; ay /= l;
const s = conv.Scale ?? 1;
return [ax * s, ay * s, 0, 0, -ay * s, ax * s, 0, 0, 0, 0, s, 0, conv.Eastings, conv.Northings, conv.OrthogonalHeight, 1];
}The angle of the model's X axis from grid east, in degrees.
θ = atan2(ordinate, abscissa) × 180 ÷ π
One atan2 divided by π ÷ 180, in 64-bit; valid in all four quadrants.
export const rotationDeg = (conv) => Math.atan2(conv.XAxisOrdinate ?? 0, conv.XAxisAbscissa ?? 1) / D2R;The direction of true north in the model's own X and Y.
g = grid direction of a 0.0001° step north; true north = (cos(g − m), sin(g − m)), m the model's rotation
Grid north is found by projecting the origin and a point 0.0001° of latitude further north, about 11 m, and taking atan2 of the difference of two grid coordinates that close together. The result is rounded to twelve decimals.
/** Where true north points in the model's own X/Y (unit vector), for IfcGeometricRepresentationContext.TrueNorth. */
export function trueNorth(conv, crs) {
try {
const ll = lonLatOf(crs, conv.Eastings, conv.Northings);
const a = gridOf(crs, ll[0], ll[1]), b = gridOf(crs, ll[0], ll[1] + 1e-4);
const g = Math.atan2(b[1] - a[1], b[0] - a[0]), m = Math.atan2(conv.XAxisOrdinate, conv.XAxisAbscissa);
return [r12(Math.cos(g - m)), r12(Math.sin(g - m))];
} catch { return null; }
}The same place on the earth described in a different coordinate system: new origin, rotation and scale.
origin' = T(E, N); p₁ = T(origin + 10 m along the X axis); axis' = (p₁ − origin') ÷ D; scale' = scale × D ÷ step
The origin goes through proj4 once. Rotation and scale come from the difference of two transformed points 10 m apart; with map coordinates in the millions that subtraction leaves about 10 of the 16 digits. Eastings, northings and height are rounded to six decimals and the axis to twelve; the scale is not rounded.
/**
* Re-express an IfcMapConversion in another projected system. The place the model stands on the earth stays the same (its latitude/longitude do
* not change); the numbers that describe it change: origin through proj4, rotation from how the new grid's north turns at the origin, scale from
* how lengths along the model's X axis change. Heights change by the unit only (a vertical datum is not touched). Nothing inside the model moves.
*/
export function reprojectConversion(conv, from, to) {
if (!from || !to || from.unit === 'deg' || to.unit === 'deg') return { error: fail('projected'), problems: [] };
const problems = [];
let before, o2;
try {
before = lonLatOf(from, conv.Eastings, conv.Northings);
o2 = transformer(from, to)(conv.Eastings, conv.Northings);
} catch { return { error: fail('outside', { name: from.name }), problems }; }
if (!before.every(fin) || !o2.every(fin)) return { error: fail('outside', { name: from.name }), problems };
problems.push(...shiftProblems(from, to), ...areaProblems(from, before[0], before[1], 'file'));
const after = lonLatOf(to, o2[0], o2[1]);
problems.push(...areaProblems(to, after[0], after[1], 'model'));
let ax = conv.XAxisAbscissa ?? 1, ay = conv.XAxisOrdinate ?? 0;
const l = Math.hypot(ax, ay) || 1; ax /= l; ay /= l;
const step = 10 / (METRES_PER[from.unit] || 1); // 10 m along the model's X axis, in the old system's unit
const p1 = transformer(from, to)(conv.Eastings + ax * step, conv.Northings + ay * step);
const dx = p1[0] - o2[0], dy = p1[1] - o2[1], D = Math.hypot(dx, dy);
if (!fin(D) || D === 0) return { error: fail('outside', { name: from.name }), problems };
const out = {
Eastings: r6(o2[0]), Northings: r6(o2[1]),
OrthogonalHeight: r6((conv.OrthogonalHeight || 0) * (METRES_PER[from.unit] || 1) / (METRES_PER[to.unit] || 1)),
XAxisAbscissa: r12(dx / D), XAxisOrdinate: r12(dy / D),
Scale: (conv.Scale ?? 1) * D / step,
};
return { conv: out, before, after, problems };
}The arithmetic is 64-bit, but the method or the data sets the accuracy: a solver tolerance, a mesh, a step size, a datum shift, a table. The limit is stated where the result is shown.