layer/geotiff/warp.tstypescript
import type { TypedArray } from "geotiff";
import { getFastNoDataChecker } from "./ColorScale";
/**
* Tile resampling data, the window rasters and the grid of the source positions of the tile pixels.
*/
export interface IWarpParams {
/** Window rasters, one by band. */
rasters: TypedArray[];
/** Source positions of the grid nodes, normalized to 0..1, v grows down the raster rows. */
gridU: Float64Array;
gridV: Float64Array;
gridSize: number;
imgWidth: number;
imgHeight: number;
window: [number, number, number, number];
tileSize: number;
bilinear: boolean;
nodata: number | null;
/** Valid pixels of the window, if the file has a mask. */
maskRaster?: TypedArray | null;
}
/**
* Returns the tile rasters resampled from the window rasters, the nearest or bilinear. Tile pixels out of
* the image, of no data or out of the mask are NaN.
*/
export function warpRasters(p: IWarpParams): Float32Array[] {
const { rasters, gridU, gridV, gridSize, imgWidth, imgHeight, tileSize, bilinear, maskRaster } = p;
const [winX0, winY0, winX1, winY1] = p.window;
const winWidth = winX1 - winX0;
const winHeight = winY1 - winY0;
const side = gridSize + 1;
const cellSize = tileSize / gridSize;
const checkNoData = p.nodata !== null ? getFastNoDataChecker(p.nodata) : null;
const out: Float32Array[] = [];
for (let b = 0; b < rasters.length; b++) {
out.push(new Float32Array(tileSize * tileSize).fill(NaN));
}
for (let y = 0; y < tileSize; y++) {
// Grid nodes sit on tile pixel borders, so pixel centers are used here
const gy = (y + 0.5) / cellSize;
const j0 = Math.min(gridSize - 1, gy | 0);
const fy = gy - j0;
for (let x = 0; x < tileSize; x++) {
const gx = (x + 0.5) / cellSize;
const i0 = Math.min(gridSize - 1, gx | 0);
const fx = gx - i0;
// Bilinear interpolation of the source position between the grid nodes
const n00 = j0 * side + i0;
const n10 = n00 + 1;
const n01 = n00 + side;
const n11 = n01 + 1;
const w00 = (1 - fx) * (1 - fy);
const w10 = fx * (1 - fy);
const w01 = (1 - fx) * fy;
const w11 = fx * fy;
const u = gridU[n00] * w00 + gridU[n10] * w10 + gridU[n01] * w01 + gridU[n11] * w11;
const v = gridV[n00] * w00 + gridV[n10] * w10 + gridV[n01] * w01 + gridV[n11] * w11;
if (u < 0 || u >= 1 || v < 0 || v >= 1) continue;
// Source position inside the window, in pixel centers
const sx = u * imgWidth - winX0 - 0.5;
const sy = v * imgHeight - winY0 - 0.5;
const outIdx = y * tileSize + x;
// Pixels the mask marks as invalid stay NaN, which is rendered as transparent
if (maskRaster) {
const mx = Math.max(0, Math.min(winWidth - 1, Math.round(sx)));
const my = Math.max(0, Math.min(winHeight - 1, Math.round(sy)));
if (!maskRaster[my * winWidth + mx]) continue;
}
if (bilinear) {
const x0 = Math.floor(sx);
const y0 = Math.floor(sy);
const tx = sx - x0;
const ty = sy - y0;
const x1 = Math.min(winWidth - 1, x0 + 1);
const y1 = Math.min(winHeight - 1, y0 + 1);
const cx0 = Math.max(0, Math.min(winWidth - 1, x0));
const cy0 = Math.max(0, Math.min(winHeight - 1, y0));
const i00 = cy0 * winWidth + cx0;
const i10 = cy0 * winWidth + x1;
const i01 = y1 * winWidth + cx0;
const i11 = y1 * winWidth + x1;
for (let b = 0; b < rasters.length; b++) {
const src = rasters[b];
const v00 = src[i00],
v10 = src[i10],
v01 = src[i01],
v11 = src[i11];
if (
checkNoData
? checkNoData(v00) || checkNoData(v10) || checkNoData(v01) || checkNoData(v11)
: Number.isNaN(v00) || Number.isNaN(v10) || Number.isNaN(v01) || Number.isNaN(v11)
) {
continue;
}
out[b][outIdx] = (v00 * (1 - tx) + v10 * tx) * (1 - ty) + (v01 * (1 - tx) + v11 * tx) * ty;
}
} else {
const px = Math.max(0, Math.min(winWidth - 1, Math.round(sx)));
const py = Math.max(0, Math.min(winHeight - 1, Math.round(sy)));
const srcIdx = py * winWidth + px;
for (let b = 0; b < rasters.length; b++) {
const val = rasters[b][srcIdx];
if (checkNoData ? checkNoData(val) : Number.isNaN(val)) continue;
out[b][outIdx] = val;
}
}
}
}
return out;
}