developmentseed / developmentseed/deck.gl-raster
Adaptive reprojection mesh never converges near the poles for global EPSG:4326 sources → high-latitude shear
- Dominant language
- TypeScript
- Stars
- 228
- Forks
- 29
- Avg merge
- 12h 27m
- Merged PRs (30d)
- 4
Description
## Summary
When a global **EPSG:4326** raster is rendered on a Web-Mercator map, tiles at high latitudes are visibly **sheared / streaked**; equatorial and mid-latitude tiles are correct. The cause is `@developmentseed/raster-reproject`'s `RasterReprojector` (used by `RasterLayer._generateMesh`): its adaptive mesh **fails to converge near the poles**, hits the `maxIterations` cap (10000), and breaks — leaving a badly under-resolved mesh for that tile.
This is already flagged in the code itself — `RasterReprojector.run()`:
```js
// Note: this primarily happens near the poles, where we'll essentially
// never converge
while (this.getMaxError() > maxError) {
this.refine();
if (++iterations > maxIterations) {
console.warn(`RasterReprojector: mesh refinement did not converge ...`);
break;
}
}
```
## Minimal repro (headless, uses this repo's own classes)
```js
// deps: @developmentseed/deck.gl-raster, @developmentseed/raster-reproject,
// @developmentseed/proj, proj4
import proj4 from "proj4";
import { AffineTilesetLevel } from "@developmentseed/deck.gl-raster";
import { RasterReprojector } from "@developmentseed/raster-reproject";
import { makeClampedForwardTo3857 } from "@developmentseed/proj";
// A plain global EPSG:4326 grid: 0.1°/px, 3600×1800, 256² tiles (7 tile rows).
const level = new AffineTilesetLevel({
affine: [0.1, 0, -180, 0, -0.1, 90], // [scaleX,0,originX,0,scaleY,originY]
arrayWidth: 3600, arrayHeight: 1800, tileWidth: 256, tileHeight: 256,
mpu: 111319.49,
});
const conv = proj4("EPSG:4326", "EPSG:3857");
const to3857 = makeClampedForwardTo3857((x, y) => conv.forward([x, y]), (x, y) => [x, y]);
const from3857 = (x, y) => conv.inverse([x, y]);
function stats(row, label) {
const { forwardTransform, inverseTransform } = level.tileTransform(0, row);
let warned = false; const w = console.warn; console.warn = () => (warned = true);
const rp = new RasterReprojector(
{ forwardTransform, inverseTransform, forwardReproject: to3857, inverseReproject: from3857 },
257, 257,
);
rp.run(0.125); // default maxError
console.warn = w;
console.log(`${label}: triangles=${rp.triangles.length / 3} finalMaxError=${rp.getMaxError().toFixed(3)} ${warned ? "NON-CONVERGENT (hit 10k cap)" : "converged"}`);
}
stats(3, "equator (-12..13°)");
stats(1, "mid-lat (39..64°) ");
stats(0, "polar (64..90°) ");
```
Output:
```
equator (-12..13°): triangles=2295 finalMaxError=0.119 converged
mid-lat (39..64°) : triangles=374 finalMaxError=0.125 converged
polar (64..90°) : triangles=10003 finalMaxError=49.489 NON-CONVERGENT (hit 10k cap)
```
The polar tile (which spans up to 90°) ends **396× over** the target error, having exhausted the iteration budget with 10k triangles that still don't fit.
## Why it never converges
`_findReprojectionCandidate` measures error in **input-pixel space**: it takes a linearly-interpolated output (Mercator) point, runs it back through `inverseReproject` (Mercator→lon/lat) and `inverseTransform` (lon/lat→pixel), and compares to the exact pixel. Near the Mercator singularity, `y(φ) = R·ln(tan(π/4 + φ/2))` diverges, so as the tile approaches 90° a bounded output error maps to an unbounded input-pixel error — the metric can't be driven below `maxError` no matter how many vertices are added. `makeClampedForwardTo3857` clamps the *forward* projection at ±85.05° (which handles the exact-pole `NaN`), but the delatin **error metric** still sees the divergent region and never settles.
Consequence: lowering `maxError` doesn't help (the polar tile never reaches any target — it always iteration-caps), which matches what we observed in a downstream app.
## Impact
Any global EPSG:4326 source drawn on a Web-Mercator viewport shears at high latitudes. It's most visible on high-contrast global data; smooth fields (e.g. temperature) hide it. Observed in a production viewer ([source-cooperative/zarr-viewer#70](https://github.com/source-cooperative/zarr-viewer/issues/70)) on a global NDVI store — equator/tropics pixel-accurate, poles sheared.
## Possible directions
- Clamp each tile's **latitude extent** to the Mercator-valid range before meshing (rows poleward of ±85.05° carry no displayable information anyway).
- Make the error metric robust to the Mercator singularity (measure in output/clip space, or cap the per-vertex Jacobian) so refinement converges.
- When the iteration cap is hit, fall back to a dense uniform grid for that tile rather than keeping the last (arbitrary) triangulation.
Happy to help test a fix against the real reproducing dataset.
Contributor guide
No contributing guide indexed for this repository
Research direction
Start with RasterLayer._generateMesh and RasterReprojector.run(), then inspect _findReprojectionCandidate and reproduce the polar case with the provided headless script. Compare latitude clamping, a robust output-space error metric, and the iteration-cap fallback before choosing an approach. Done means global EPSG:4326 tiles near the poles converge or use a defined fallback without shear, warning, or excessive residual error.
Written by the indexing model from the issue text.
Assessment
- Tech stack
- typescript
- Domain
- computer-graphics, frontend
- Issue type
- Bug
- Difficulty
- 5/5
- Estimated time
- Over a week
- Activity status
- Quiet
- Clarity
- Mostly clear
- Newbie friendliness
- 42/100