PR #349: Antimeridian longitude wrap + NaN from oversized tile matrices on global COGs

[yharby]

Context

PR #349 (merged as 6adddc4) introduces CARTESIAN/EPSG:3857 rendering for WebMercator viewports. This works great for regional COGs but produces rendering artifacts on large global EPSG:4326 datasets where tiles span ±180° longitude.

Test file

gdalinfo /vsicurl/https://s3.us-west-2.amazonaws.com/us-west-2.opendata.source.coop/walkthru-earth/dem-terrain/GEDTM30/gedtm30.tif

1440010×600010 pixels, EPSG:4326, Int32, Deflate, 2048×2048 tiles, 10 overviews, bbox [-180.00125, -65.00125, 180.00125, 85.00125].

Issues found

1. inverseFrom3857 returns +180° at the antimeridian → 360° round-trip error

For tiles touching lon=-180°, the mesh refinement round-trip CRS → 3857 → CRS can flip longitude from -180° to +180° due to floating-point ambiguity in proj4. This causes a round-trip error equal to the full image width (5625 pixels for this file), and the mesh never converges.

Workaround: Normalize longitude in inverseFrom3857/inverseFrom4326:

if (r[0] >= 180) r[0] -= 360;
else if (r[0] < -180) r[0] += 360;

2. sampleReferencePointsInEPSG3857 receives NaN from oversized tile matrices

When tileHeight > imageHeight in a COG overview (e.g. 2048px tile for a 2343px overview), the tile matrix bounds extend far beyond the data bbox (lat=-177° for this file). projectTo3857 returns NaN for these out-of-domain coordinates, which propagates through tile traversal → "Invalid number null" crash in TileLayer initialization.

Workaround: Clamp non-finite values in sampleReferencePointsInEPSG3857:

if (!Number.isFinite(projected[0])) projected[0] = 0;
if (!Number.isFinite(projected[1])) projected[1] = geoY > 0 
  ? EPSG_3857_HALF_CIRCUMFERENCE : -EPSG_3857_HALF_CIRCUMFERENCE;

3. RasterReprojector.run() has no iteration limit → browser hang

For tiles with extreme Mercator non-linearity (partial tiles spanning -65° to -177° latitude), the Delatin mesh refinement while (getMaxError() > maxError) { refine() } loop never converges and hangs the browser.

Workaround: Add iteration cap in delatin.js:

const MAX_ITERATIONS = 5000;
while (this.getMaxError() > maxError) {
  this.refine();
  if (++iterations >= MAX_ITERATIONS) break;
}

4. Partial edge tiles produce distorted meshes

At z=0, tiles where actual data height is <25% of tile matrix height (e.g. 295px data in 2048px tile) have forwardTransform mapping UV [0,1] to CRS coordinates far beyond the data bbox. The mesh is built for the full UV range, creating distorted geometry.

Workaround: Skip rendering tiles where data.height / tileMatrix.tileHeight < 0.25.

5. Bbox slightly beyond ±180° wraps easting sign

The test file has bbox lon=±180.00125°. forwardTo3857(-180.00125, lat) wraps to positive easting (+20037369 instead of -20037647), placing the tile on the wrong side of the world.

Workaround: Clamp EPSG:4326 bbox to exactly [-180, lat_min, 180, lat_max] before calling generateTileMatrixSet.

Environment

• @developmentseed/deck.gl-geotiff@0.3.0 + PR fix: Render mesh from Web Mercator coordinates #349 backport
• @developmentseed/deck.gl-raster@0.3.0
• @developmentseed/raster-reproject@0.3.0
• deck.gl 9.2, MapLibre 5.20, Firefox/Chrome on macOS

[yharby]

note that after i zoom in to preview the next level of overviews it works soo good and the projection is correct

[kylebarron]

We don't currently claim to support images spanning the antimeridian. See #171 and #269.

Regardless, switching to rendering from EPSG:3857 is crucial for correctness.

3. RasterReprojector.run() has no iteration limit → browser hang

This isn't true on latest main.