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

[tylere]

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():

// 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)

// 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) 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.

[kylebarron]

I think this should already be fixed. I thought it was already released but perhaps not (I've been focused on other things and I need to release 0.8

Oh, indeed, this is fixed on main but not released
https://github.com/developmentseed/deck.gl-raster/pull/618/changes#diff-61a40941a8bfa8b1b6e35c973d994cf0870877717170de55ef13ec4b2c433d87R50