Reads a terrain grid's position, cell size and elevations from a GeoTIFF or an ESRI ASCII grid, with its scale, offset and no-data value.
How the precision is kept
The grid's position comes from the file's georeferencing tags in 64-bit numbers, and elevations are stored in a Float64Array after scale and offset are applied. The digits an elevation has are those of the file: a 32-bit float GeoTIFF carries about 7, an ASCII grid as many as were written. Reading at a step keeps every step-th cell without averaging.
Calculations under this policy
8 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.
A ModelTransformation matrix is copied as it is; otherwise the origin is one multiplication and one subtraction or addition from the first tie point and the pixel scale. A file marked pixel-is-point is shifted by half a cell so the origin is a cell corner. A transform whose determinant is 0 is refused.
raster.js
lines 105 to 112
BimCraft 0.8.1
if (mtx && mtx.length >= 16 && mtx.every(Number.isFinite)) gt = [mtx[3], mtx[0], mtx[1], mtx[7], mtx[4], mtx[5]];
else if (tie && tie.length >= 6 && scl && scl.length >= 2 && [...tie.slice(0, 6), ...scl.slice(0, 2)].every(Number.isFinite) && scl[0] > 0 && scl[1] > 0) {
if (tie.length > 6) notes.push(t('raster.tie_points'));
gt = [tie[3] - tie[0] * scl[0], scl[0], 0, tie[4] + tie[1] * scl[1], 0, -scl[1]];
} else throw new Error(t('raster.no_georef'));
if (!(Math.abs(gt[1] * gt[5] - gt[2] * gt[4]) > 0)) throw new Error(t('raster.zero_cell'));
// pixel-is-point: the tie point is a cell center (GDAL's reading) — shift to the corner convention
if (pixelIsPoint) { gt = [gt[0] - 0.5 * (gt[1] + gt[2]), gt[1], gt[2], gt[3] - 0.5 * (gt[4] + gt[5]), gt[4], gt[5]]; }
Cell size
64-bit
The single cell size of a square, north-up grid.
cell = a when b = d = 0 and |a + e| < 1e-9 × |a|
The squareness test is relative to the cell size. A rotated or non-square grid has no single cell size and the value is left empty.
The size, cell size and lower-left corner of an ASCII grid.
corner = center − cell ÷ 2 when the header gives cell centers
Header numbers are parsed with Number. A header that gives xllcenter is shifted by half a cell in x and y; if it gives no yllcenter, 0 is used for it. Column and row counts must be whole numbers.
raster.js
lines 135 to 155
BimCraft 0.8.1
async function ascHeader(file) {
const head = await file.slice(0, 4096).text();
const lines = head.split(/\r?\n/);
const h = {};
let bytes = 0, n = 0;
for (const line of lines) {
const m = line.match(/^\s*([A-Za-z_]+)\s+(\S+)\s*$/);
if (!m || !ASC_KEYS.includes(m[1].toLowerCase())) break;
h[m[1].toLowerCase()] = Number(m[2]);
bytes += new TextEncoder().encode(line).length + (head.includes('\r\n') ? 2 : 1);
if (++n > 10) break;
}
const ncols = h.ncols, nrows = h.nrows;
if (!Number.isInteger(ncols) || !Number.isInteger(nrows) || ncols <= 0 || nrows <= 0 || ncols > LIMITS.dim || nrows > LIMITS.dim) throw new Error(t('raster.asc.not_asc'));
const dx = h.cellsize ?? h.dx, dy = h.cellsize ?? h.dy ?? h.dx;
if (!(dx > 0) || !(dy > 0) || !Number.isFinite(dx) || !Number.isFinite(dy)) throw new Error(t('raster.asc.no_cell'));
const center = h.xllcenter !== undefined;
const xll = center ? h.xllcenter - dx / 2 : h.xllcorner, yll = center ? (h.yllcenter ?? 0) - dy / 2 : h.yllcorner;
if (!Number.isFinite(xll) || !Number.isFinite(yll)) throw new Error(t('raster.asc.no_corner'));
return { ncols, nrows, dx, dy, xll, yll, nodata: Number.isFinite(h.nodata_value) ? h.nodata_value : null, dataStart: bytes };
}
ASCII grid position
64-bit
The transform of an ASCII grid, with its origin at the top-left corner.
top = yll + rows × dy
The top edge is the lower-left y plus the row count times the cell height, one multiplication and one addition. The grid has a single cell size only when dx equals dy exactly.
Turns one stored sample into an elevation, or into no value.
z = v × scale + offset
The multiplication and addition are 64-bit. No-data is matched by exact equality; for a 32-bit float file the sample is also compared with Math.fround of the no-data value, because such files store the 32-bit rounding of it. The sample itself has only the digits of the file's type.
raster.js
lines 207 to 213
BimCraft 0.8.1
function cleanValue(src, v) {
if (src.nodata != null && (v === src.nodata || (Number.isNaN(src.nodata) && Number.isNaN(v)))) return NaN;
if (!Number.isFinite(v)) return NaN;
// float32 files often store nodata as the float32 rounding of a double (−3.4028235e38): compare at float32 too
if (src.nodata != null && src.dtype === 'float32' && Math.fround(src.nodata) === v) return NaN;
return v * src.scale + src.offset;
}
Reading a window of cells
64-bit
Reads part of the grid, optionally every step-th cell, into an array of elevations.
The window is widened to whole cells with Math.floor and Math.ceil and the step is truncated to a whole number of at least 1. Cells between steps are skipped, not averaged. Values are kept in a Float64Array.
The minimum and maximum elevation shown for a grid.
Statistics stored in the file are used when present, with scale and offset applied. Otherwise a GeoTIFF is read at 256 × 256 by nearest neighbor, so a grid larger than that can have extremes the sample misses; the result is marked as sampled. An ASCII grid up to 20 MB is read in full.
raster.js
lines 281 to 297
BimCraft 0.8.1
export async function sampleStats(src, { signal } = {}) {
if (src.stats) return src.stats;
let min = Infinity, max = -Infinity, sampled = true;
if (src.kind === 'geotiff') {
const w = Math.min(256, src.width), h = Math.min(256, src.height);
const [data] = await src.image.readRasters({ samples: [0], width: w, height: h, resampleMethod: 'nearest', signal: signal || undefined });
for (let i = 0; i < data.length; i++) { const v = cleanValue(src, data[i]); if (v < min) min = v; if (v > max) max = v; }
sampled = w < src.width || h < src.height;
} else {
if (src.size > 20e6) return null;
const g = await readWindow(src, { signal });
for (const v of g.z) { if (v < min) min = v; if (v > max) max = v; }
sampled = false;
}
src.stats = Number.isFinite(min) ? { min, max, sampled } : { min: null, max: null, sampled };
return src.stats;
}
Position of a georeferenced picture
64-bit
The transform and bounding box of a picture read at a reduced size.
When a picture is reduced, its pixel counts are truncated with Math.floor and the cell terms are scaled by the exact ratio of the original to the reduced counts, so the picture's outer edges stay where they were. The box is the minimum and maximum of the four transformed corners.
Affine transform from ModelTransformation, or from the first tie point and the pixel scale; half-cell shift for pixel-is-point files; value × scale + offset; decimation by a whole step. GeoTIFF decoding is done by the geotiff library the file imports.
Number types
64-bit numbers and Float64Array for what is kept; the file's own sample type (8, 16 or 32-bit integer, 32 or 64-bit float, or text) for what is read.
Tolerances
A grid counts as square and north-up when |x cell + y cell| < 1e-9 × |x cell| (relative) with no rotation terms. No-data is matched by exact equality, and for 32-bit float files also against the 32-bit rounding of the no-data value.
Limits
Only the first tie point is used when a file has several (a note says so). The minimum and maximum of a GeoTIFF without stored statistics are taken from a 256 × 256 nearest-neighbor sample and marked as sampled.
Approach: 64-bit, stated limit
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.