developmentseed / developmentseed/deck.gl-raster

Adaptive reprojection mesh never converges near the poles for global EPSG:4326 sources → high-latitude shear

Open
#625 1 comment 0 reactions 0 assignees View on GitHub
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

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.