Resources Numerical accuracy Calculation policy
Converts project coordinates to and from a named map projection, longitude and latitude, and other systems, and reports grid scale factor, convergence and elevation factor at a point.
Project coordinates stay 64-bit numbers. The optional local step (scale, rotation, two base points) is a closed-form similarity. The projection and datum shift are done by the proj4 library the file imports, so their accuracy is that library's and that of the datum shift in the definition. Scale factor and convergence are finite differences with a fixed step of 1e-5 radians. A definition is matched to a catalog system when three probe points agree to 1e-3 of the system's unit.
17 calculations in 2 source files. 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 bin/test-civilcraft-crs.mjs, bin/test-crs-share.mjs
The length of a meter, an international foot and a US survey foot in meters.
m = 1; ft = 0.3048; ftUS = 1200 ÷ 3937
0.3048 is the defined value of the international foot. The US survey foot is written as its defining fraction 1200 ÷ 3937 and divided once in 64-bit, not typed as a rounded decimal.
export const METRES_PER = { m: 1, ft: 0.3048, ftUS: 1200 / 3937 };Moves a point between project (ground or assumed) coordinates and grid coordinates.
grid = to + k × R(θ) × (p − from); p = from + R(θ)ᵀ × (grid − to) ÷ k
The rotation is given in degrees and converted with Math.PI ÷ 180; cosine and sine are taken once. The offset from the base point is subtracted first, then rotated and scaled, so the terms combined are small differences and not full coordinates. The inverse divides by k without a check; crsProblem refuses a scale that is not above zero.
const c = Math.cos(l.rotation * Math.PI / 180), s = Math.sin(l.rotation * Math.PI / 180), k = l.scale;
toGrid = (x, y) => { const dx = x - l.from[0], dy = y - l.from[1]; return [l.to[0] + k * (c * dx - s * dy), l.to[1] + k * (s * dx + c * dy)]; };
fromGrid = (gx, gy) => { const dx = (gx - l.to[0]) / k, dy = (gy - l.to[1]) / k; return [l.from[0] + c * dx + s * dy, l.from[1] - s * dx + c * dy]; };The longitude and latitude (WGS 84) of a project point, and the project point of a longitude and latitude.
The local step runs first, then proj4 converts between the system's definition and WGS 84. With a loaded grid file the definition's shift is replaced by the grid (defFor); otherwise the shift written in the definition is used. This code adds no rounding of its own.
toLonLat(x, y) { const g = toGrid(x, y); return geo.forward(g); },
fromLonLat(lon, lat) { const g = geo.inverse([lon, lat]); return fromGrid(g[0], g[1]); },At a project point: the grid scale factor, the meridian convergence in degrees, the elevation factor and the combined factor.
k_N = |p(φ+d) − p(φ−d)| ÷ (2 M d); k_E = |p(λ+d) − p(λ−d)| ÷ (2 N cos φ d); k = (k_N + k_E) ÷ 2; γ = −atan2(Δx, Δy); elevation = R ÷ (R + h), R = √(M N)
Central differences with a fixed step d = 1e-5 radians: the result carries a truncation error of order d² and the rounding of a difference of two projected coordinates, so it does not hold 15 digits. M and N use the system's own ellipsoid. h is z converted to meters and taken as ellipsoid height. The division by cos φ is not guarded at the poles.
/** At a project point: grid scale factor k (along the meridian and the parallel, averaged — conformal systems make
* them equal), convergence γ (degrees, the surveying sign: geodetic azimuth = grid azimuth + γ; + east of the
* central meridian in the northern hemisphere; PROJ's meridian_convergence),
* elevation factor R / (R + h) with h = z taken as ellipsoid height (no geoid model: off by ~5 ppm per 30 m of
* geoid separation — say so where it is shown), combined = k × elevation. All on the system's own ellipsoid. */
factors(x, y, z = 0) {
const P = own(), toM = P.to_meter || 1;
const g = toGrid(x, y);
const ll = P.inverse({ x: g[0] * toM, y: g[1] * toM });
const lam = ll.x, phi = ll.y, d = 1e-5; // radians (~60 m; central differences: the error is ~d², rounding ~1e-10)
// central differences (builder C, crs tests: one-sided ones were off by up to 2e-7 in k and 4e-6° in γ vs PROJ)
const pN = P.forward({ x: lam, y: phi + d }), pS = P.forward({ x: lam, y: phi - d }), pE = P.forward({ x: lam + d, y: phi }), pW = P.forward({ x: lam - d, y: phi });
const e2 = P.es || 0, w = Math.sqrt(1 - e2 * Math.sin(phi) ** 2);
const N = P.a / w, M = P.a * (1 - e2) / w ** 3;
const kN = Math.hypot(pN.x - pS.x, pN.y - pS.y) / (2 * M * d), kE = Math.hypot(pE.x - pW.x, pE.y - pW.y) / (2 * N * Math.cos(phi) * d);
const k = (kN + kE) / 2;
const convergence = -Math.atan2(pN.x - pS.x, pN.y - pS.y) * 180 / Math.PI;
const R = Math.sqrt(M * N);
const h = (Number(z) || 0) * (METRES_PER[crs.unit] || 1);
const elevation = R / (R + h);
const [lon, lat] = geo.forward(g);
return { k, kMeridian: kN, kParallel: kE, convergence, elevation, combined: k * elevation, lon, lat };
},The scale and rotation that carry two points known in both local and grid coordinates onto each other.
scale = |B − A| ÷ |b − a|; rotation = atan2(grid) − atan2(local), in degrees, brought into −180 to 180
Lengths by Math.hypot and angles by Math.atan2 of coordinate differences. Both spans must be longer than 1e-6 (absolute, in the coordinates' own unit). With exactly two points there is no redundancy, so an error in either point goes straight into the scale and rotation.
/** Two control points [[local, grid], [local, grid]] (each [x, y]) → `local` (scale, rotation about the first). */
export function localFromPairs(pairs) {
if (!Array.isArray(pairs) || pairs.length !== 2) throw new Error(t('crs.pairs.two'));
const [[a, A], [b, B]] = pairs;
const dl = [b[0] - a[0], b[1] - a[1]], dg = [B[0] - A[0], B[1] - A[1]];
const ll = Math.hypot(...dl), lg = Math.hypot(...dg);
if (!(ll > 1e-6) || !(lg > 1e-6)) throw new Error(t('crs.pairs.same'));
const rotation = (Math.atan2(dg[1], dg[0]) - Math.atan2(dl[1], dl[0])) * 180 / Math.PI;
return { scale: lg / ll, rotation: ((rotation + 540) % 360) - 180, from: [a[0], a[1]], to: [A[0], A[1]] };
}A local system that scales by a combined factor about one base point, with no rotation.
scale = combined; rotation = 0; from = to = base
The factor and the base point are stored as given; no arithmetic is done here.
/** Ground coordinates: scale about a base point (the same point in both systems). */
export const groundLocal = (combined, base) => ({ scale: combined, rotation: 0, from: [base[0], base[1]], to: [base[0], base[1]] });Converts a coordinate pair from one system's definition to another's.
One proj4 forward call between the two definitions; a missing system means WGS 84 longitude and latitude. The local step is not applied here. Converters are cached by their two definitions.
/** Any two systems (state.crs shapes, catalog entries, proj4/WKT strings; null = WGS 84 lon/lat). Ignores `local`. */
export function transformer(from, to) {
const t = converter(defOf(from), defOf(to));
return (x, y) => t.forward([x, y]);
}Decides whether the projection in a .prj file is one of the catalog's systems, and if so uses the catalog's definition.
|a − a′| < 1e-3; |e² − e′²| < 1e-9; distance at each of 3 probe points < 1e-3 (1e-9 for geographic)
Both definitions are projected at three probe points on their own ellipsoids and compared with absolute tolerances. A match replaces the file's definition with the catalog's, which may carry a datum shift the file left out. The 1e-3 tolerance is in the system's unit, so it is 1 mm for meters and about 0.3 mm for feet.
const sameEllipsoid = (P) => Math.abs(P.a - def.a) < 1e-3 && Math.abs((P.es || 0) - (def.es || 0)) < 1e-9; // a .prj rounds 1/f
// geographic: radians either way (a WKT's degree UNIT sets to_meter ≈ 111 km — builder C: GEOGCS .prj files never matched)
const fwd = (P, [lon, lat]) => { const r = P.forward({ x: lon * Math.PI / 180, y: lat * Math.PI / 180 }), k = geographic ? 1 : (P.to_meter || 1); return [r.x / k, r.y / k]; };
let want = null;
try { want = pts.map((p) => fwd(def, p)); } catch { want = null; }
if (want && want.every((w) => w.every(Number.isFinite))) {
const tol = geographic ? 1e-9 : 1e-3;
for (const e of CATALOG) {
if (e.unit !== unit || e.geographic !== geographic) continue;
try {
const P = proj4.Proj(e.proj4);
if (!sameEllipsoid(P)) continue;
if (pts.every((p, i) => { const g = fwd(P, p); return Math.hypot(g[0] - want[i][0], g[1] - want[i][1]) < tol; })) {
if (!matched || prjRank(e, text) < prjRank(matched, text)) matched = e;
}Three longitude and latitude points inside a projection's zone, used to compare two definitions.
(λ₀, φ), (λ₀ + 0.5°, φ + 0.5°), (λ₀ − 0.5°, φ + 0.25°)
The points lie on and half a degree off the central meridian, at the mean of the two standard parallels or the latitude of origin, with the latitude kept within ±80 degrees. Two definitions that agree at these three points but differ elsewhere would be taken as the same.
// three points inside the zone, as lon/lat (for matching a .prj to the catalog): on the central meridian, between the
// standard parallels when there are two, else at the latitude of origin (UTM: 45°, so the false northing is exercised)
function probePoints(def, geographic) {
const deg = 180 / Math.PI;
const lon0 = def.zone ? def.zone * 6 - 183 : (def.long0 || 0) * deg;
let lat = def.lat1 != null && def.lat2 != null ? (def.lat1 + def.lat2) / 2 * deg : (def.lat0 || 0) * deg;
if (def.zone) lat = def.utmSouth ? -45 : 45;
if (geographic) lat = 0;
lat = Math.max(-80, Math.min(80, lat));
return [[lon0, lat], [lon0 + 0.5, lat + 0.5], [lon0 - 0.5, lat + 0.25]];
}Decides whether a projection is in meters, international feet or US survey feet.
|k − 1| < 1e-12 → m; |k − 0.3048| < 1e-10 → ft; |k − 1200 ÷ 3937| < 1e-10 → ftUS
proj4's unit name is used when it gives one. Otherwise the factor in the UNIT clause (or proj4's to_meter) is compared with absolute tolerances, so a US survey foot factor rounded to 10 decimals is still recognized. Any other factor is refused.
function unitFromDef(def, text) {
const u = def.units && UNIT_OF_PROJ4[def.units];
if (u) return u;
const m = text.match(/UNIT\["([^"]+)",\s*([\d.eE+-]+)\]\s*(,\s*AXIS|,\s*AUTHORITY|\]\s*$)/i) || text.match(/UNIT\["([^"]+)",\s*([\d.eE+-]+)/gi)?.slice(-1)[0]?.match(/UNIT\["([^"]+)",\s*([\d.eE+-]+)/i);
const k = m ? Number(m[2]) : (def.to_meter || 1);
if (Math.abs(k - 1) < 1e-12) return 'm';
if (Math.abs(k - 0.3048) < 1e-10) return 'ft';
if (Math.abs(k - 1200 / 3937) < 1e-10) return 'ftUS';
return null;
}The west, south, east and north limits, in degrees, of the area a loaded grid shift file covers.
east = −λ_ll × 180 ÷ π; west = −(λ_ll + (n_λ − 1) × Δλ) × 180 ÷ π; north = (φ_ll + (n_φ − 1) × Δφ) × 180 ÷ π
Computed from each subgrid's corner, cell size and cell count as proj4 holds them (radians, longitude positive west), then the minimum and maximum over all subgrids. The extent is used to say whether a point lies inside the grid.
const D = 180 / Math.PI;
let w = Infinity, s = Infinity, e = -Infinity, n = -Infinity;
for (const sg of g.subgrids) {
const east = -sg.ll[0] * D, west = -(sg.ll[0] + (sg.lim[0] - 1) * sg.del[0]) * D, south = sg.ll[1] * D, north = (sg.ll[1] + (sg.lim[1] - 1) * sg.del[1]) * D;
w = Math.min(w, west); e = Math.max(e, east); s = Math.min(s, south); n = Math.max(n, north);
}Tests bin/test-crs-share.mjs
The EPSG number of a UTM zone typed as text such as UTM 33N.
WGS 84: 32600 + zone (north), 32700 + zone (south); NAD83: 26900 + zone; ETRS89: 25800 + zone; GDA94: 28300 + zone; GDA2020: 7800 + zone
Whole-number arithmetic on a zone from 1 to 60. A zone outside the range the code lists for a datum returns null.
/** "UTM 33N" etc. → EPSG number (WGS 84 unless a datum is named), or null. */
export function utmCode(text) {
const s = String(text ?? '').trim();
const m = /^(?:(wgs\s*-?\s*84|nad\s*-?\s*83|etrs\s*-?\s*89|gda\s*-?\s*94|gda\s*-?\s*2020)\s*[/ ]?\s*)?(?:utm\s*(?:zone)?\s*)(\d{1,2})\s*(n|s|north|south)?\s*$/i.exec(s);
if (!m) return null;
const zone = Number(m[2]);
if (zone < 1 || zone > 60) return null;
const south = /^s/i.test(m[3] || '');
const d = (m[1] || 'wgs84').toLowerCase().replace(/[\s-]/g, '');
if (d === 'wgs84') return (south ? 32700 : 32600) + zone;
if (d === 'nad83') return !south && zone >= 1 && zone <= 23 ? 26900 + zone : null;
if (d === 'etrs89') return !south && zone >= 28 && zone <= 38 ? 25800 + zone : null;
if (d === 'gda94') return south && zone >= 48 && zone <= 58 ? 28300 + zone : null;
if (d === 'gda2020') return south && zone >= 46 && zone <= 59 ? 7800 + zone : null;
return null;
}Decides whether a typed projection definition is in meters, international feet, US survey feet or degrees.
|k − 1| < 1e-12 → m; |k − 0.3048| < 1e-10 → ft; |k − 1200 ÷ 3937| < 1e-10 → ftUS
A +units name is looked up; otherwise +to_meter is parsed as a number and compared with the same absolute tolerances crs.js uses. Any other factor is refused.
function unitOfProj4(def, text) {
if (/\+proj=longlat|\+proj=latlong/i.test(text)) return 'deg';
const u = /\+units=(\S+)/i.exec(text);
if (u) return UNIT_OF[u[1].toLowerCase()] || null;
const k = /\+to_meter=([\d.eE+-]+)/i.exec(text);
const m = k ? Number(k[1]) : 1;
if (Math.abs(m - 1) < 1e-12) return 'm';
if (Math.abs(m - 0.3048) < 1e-10) return 'ft';
if (Math.abs(m - 1200 / 3937) < 1e-10) return 'ftUS';
return null;
}Three longitude and latitude points inside the definition's zone, used to compare it with catalog systems.
The same three points as crs.js probePoints, without the geographic case: on the central meridian and half a degree to each side, latitude kept within ±80 degrees.
// points on the zone's central meridian and a little off it, as crs.js does for a .prj
function probes(def) {
const deg = 180 / Math.PI;
const lon0 = def.zone ? def.zone * 6 - 183 : (def.long0 || 0) * deg;
let lat = def.lat1 != null && def.lat2 != null ? (def.lat1 + def.lat2) / 2 * deg : (def.lat0 || 0) * deg;
if (def.zone) lat = def.utmSouth ? -45 : 45;
lat = Math.max(-80, Math.min(80, lat));
return [[lon0, lat], [lon0 + 0.5, lat + 0.5], [lon0 - 0.5, lat + 0.25]];
}Finds the catalog systems that give the same easting and northing as a typed definition.
same projection family; |a − a′| ≤ 1e-3; |e² − e′²| ≤ 1e-9; distance at each of 3 probe points < 1e-3
Absolute tolerances: 1e-3 of the system's unit at three points, on each definition's own ellipsoid. Every catalog system that passes is kept and ranked by name and datum; the first becomes the system used.
/** The catalog system that gives the same easting / northing as this definition (own ellipsoid, 1 mm), or null. */
function matchCatalog(proj, unit, hint = '') {
if (!catalogLoaded() || unit === 'deg') return null;
let def;
try { def = proj4.Proj(proj); } catch { return null; }
const pts = probes(def);
let want;
try { want = pts.map((p) => fwd(def, p)); } catch { return null; }
if (!want.every((w) => w.every(Number.isFinite))) return null;
const family = (text) => { const k = /\+proj=(\w+)/.exec(text)?.[1]; return k === 'utm' || k === 'etmerc' || k === 'tmerc' ? 'tm' : k; }; // UTM is a transverse Mercator
const kind = family(proj), mine = datumTokens(proj);
const all = [];
for (const e of catalogEntries()) {
if (e.unit !== unit || e.geographic) continue;
if (kind && family(e.proj4) !== kind) continue;
try {
const P = proj4.Proj(e.proj4);
if (Math.abs(P.a - def.a) > 1e-3 || Math.abs((P.es || 0) - (def.es || 0)) > 1e-9) continue;
if (pts.every((p, i) => { const g = fwd(P, p); return Math.hypot(g[0] - want[i][0], g[1] - want[i][1]) < 1e-3; })) all.push(e);
} catch { /* a definition proj4js can't build is not a match */ }
}
if (!all.length) return null;
all.sort((a, b) => rankOf(a, mine, hint) - rankOf(b, mine, hint));
return { best: all[0], others: all.slice(1) };
}Refuses a typed projection that does not return its own probe point.
|inverse(forward(λ, φ)) − (λ, φ)| < 1e-6 degrees in each of λ and φ
One point is projected and brought back through proj4; each angle must return within 1e-6 degrees (absolute), about 0.1 m on the ground. It is a usability check on the definition, not a statement of the projection's accuracy.
const [lon, lat] = probes(def)[0], fw = proj4('EPSG:4326', s), q = fw.forward([lon, lat]), b = fw.inverse(q);
if (!q.every(Number.isFinite) || !(Math.abs(b[0] - lon) < 1e-6 && Math.abs(b[1] - lat) < 1e-6)) throw new Error('roundtrip');Says whether two systems give the same easting and northing.
same unit and distance at each of 3 probe points < 1e-3
Three probe points of the first system are projected through both definitions and must agree within 1e-3 of the unit (absolute). The ellipsoids are not compared here.
export function sameCoordinates(a, b) {
if (!a || !b || a.unit !== b.unit) return false;
try {
const A = proj4.Proj(a.proj4), B = proj4.Proj(b.proj4);
const pts = probes(A);
return pts.every((p) => { const x = fwd(A, p), y = fwd(B, p); return Math.hypot(x[0] - y[0], x[1] - y[1]) < 1e-3; });
} catch { return false; }
}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.