fix: render COGs in bounded projections (Mollweide, Sinusoidal, Equal Earth)

packages/deck.gl-raster/src/raster-layer.ts

@@ -247,6 +247,68 @@ export class RasterLayer extends CompositeLayer<RasterLayerProps> {
     // To account for this, we add 1 to both width and height when generating
     // the mesh. This also solves obvious gaps in between neighboring tiles in
     // the COGLayer.
+    //
+    // For tiles that straddle the CRS domain boundary (e.g. a Mollweide tile
+    // with one corner outside the ellipse), the adaptive reprojector wastes its
+    // entire iteration budget refining OOD-vertex triangles, leaving the valid
+    // area under-refined (reprojection error 100s of pixels). Detect these
+    // "border tiles" up front and use a dense uniform grid instead: OOD-vertex
+    // triangles are simply filtered out, producing a clean domain-edge cutoff
+    // at predictable 1/(BORDER_GRID_SIZE) tile-fraction resolution.
+    const { forwardTransform, forwardReproject, inverseReproject } =
+      reprojectionFns;
+    const borderCorners: [number, number][] = [
+      [0, 0],
+      [1, 0],
+      [0, 1],
+      [1, 1],
+    ];
+    // Tile width in projected CRS units for the round-trip error threshold.
+    const tileWidthCRS = Math.abs(
+      forwardTransform(width, 0)[0] - forwardTransform(0, 0)[0],
+    );
+    const roundTripThreshold = tileWidthCRS * 0.01;
+    const isBorderTile = borderCorners.some(([u, v]) => {
+      const [ix, iy] = forwardTransform(u * width, v * height);
+      const [ox, oy] = forwardReproject(ix, iy);
+      if (!Number.isFinite(ox) || !Number.isFinite(oy)) {
+        return true;
+      }
+      // Round-trip check: OOD corners where forwardReproject clamps to a
+      // finite domain-boundary position instead of returning NaN. The clamped
+      // value re-projects back to a CRS position far from the original,
+      // revealing the OOD via large round-trip error.
+      if (inverseReproject && roundTripThreshold > 0) {
+        const [ix2, iy2] = inverseReproject(ox, oy);
+        if (!Number.isFinite(ix2) || !Number.isFinite(iy2)) {
+          return true;
+        }
+        const err = Math.sqrt(
+          (ix2 - ix) * (ix2 - ix) + (iy2 - iy) * (iy2 - iy),
+        );
+        if (err > roundTripThreshold) {
+          return true;
+        }
+      }
+      return false;
+    });
+    if (isBorderTile) {
+      const { indices, positions64High, positions64Low, texCoords } =
+        buildClippedGridMesh(reprojectionFns, width + 1, height + 1);
+      this.setState({
+        reprojector: undefined,
+        mesh: {
+          indices: { value: indices, size: 1 },
+          attributes: {
+            POSITION: { value: positions64High, size: 3 },
+            TEXCOORD_0: { value: texCoords, size: 2 },
+          },
+        },
+        positions64Low,
+      });
+      return;
+    }
+
     const reprojector = new RasterReprojector(
       reprojectionFns,
       width + 1,
@@ -379,19 +441,38 @@ function reprojectorToMesh(reprojector: RasterReprojector): {
   const texCoords = new Float32Array(reprojector.uvs);
 
   const positions = new Float64Array(numVertices * 3);
+  // Track which vertices are outside the CRS domain (NaN output position).
+  const isOOD = new Uint8Array(numVertices);
   for (let i = 0; i < numVertices; i++) {
-    positions[i * 3] = reprojector.exactOutputPositions[i * 2]!;
-    positions[i * 3 + 1] = reprojector.exactOutputPositions[i * 2 + 1]!;
+    const x = reprojector.exactOutputPositions[i * 2]!;
+    const y = reprojector.exactOutputPositions[i * 2 + 1]!;
+    positions[i * 3] = x;
+    positions[i * 3 + 1] = y;
     // z (flat on the ground)
     positions[i * 3 + 2] = 0;
+    if (!Number.isFinite(x) || !Number.isFinite(y)) {
+      isOOD[i] = 1;
+    }
+  }
+
+  // Filter out any triangle that contains an out-of-domain vertex. This
+  // clips the rendered mesh cleanly at the CRS boundary without relying on
+  // undefined GPU NaN behaviour.
+  const allTriangles = reprojector.triangles;
+  const filteredTriangles: number[] = [];
+  for (let t = 0; t * 3 < allTriangles.length; t++) {
+    const a = allTriangles[t * 3]!;
+    const b = allTriangles[t * 3 + 1]!;
+    const c = allTriangles[t * 3 + 2]!;
+    if (!isOOD[a] && !isOOD[b] && !isOOD[c]) {
+      filteredTriangles.push(a, b, c);
+    }
   }
 
   // Split the float64 positions into high and low parts for fp64 emulation in
   // the shader.
   const [positions64Low, positions64High] = splitFloat64Array(positions);
-
-  // TODO: Consider using 16-bit indices if the mesh is small enough
-  const indices = new Uint32Array(reprojector.triangles);
+  const indices = new Uint32Array(filteredTriangles);
 
   return {
     indices,
@@ -400,3 +481,96 @@ function reprojectorToMesh(reprojector: RasterReprojector): {
     texCoords,
   };
 }
+
+/**
+ * Build a dense uniform grid mesh over a tile, filtering out any triangle
+ * that contains a vertex whose output position is outside the CRS domain
+ * (i.e. forwardReproject returned NaN or fails the round-trip check). Used
+ * for "border tiles" where the CRS boundary passes through the tile; the
+ * adaptive reprojector handles these poorly because OOD-vertex triangles
+ * consume its entire iteration budget and leave the valid area under-refined.
+ */
+const BORDER_GRID_SIZE = 64;
+
+function buildClippedGridMesh(
+  reprojectionFns: ReprojectionFns,
+  width: number,
+  height: number,
+  gridSize = BORDER_GRID_SIZE,
+): {
+  indices: Uint32Array;
+  positions64High: Float32Array;
+  positions64Low: Float32Array;
+  texCoords: Float32Array;
+} {
+  const { forwardTransform, forwardReproject, inverseReproject } =
+    reprojectionFns;
+  const cols = gridSize;
+  const rows = gridSize;
+  const numVerts = (cols + 1) * (rows + 1);
+
+  const positions = new Float64Array(numVerts * 3);
+  const texCoords = new Float32Array(numVerts * 2);
+  const isOOD = new Uint8Array(numVerts);
+
+  // Tile CRS width for round-trip error threshold (same logic as isBorderTile).
+  const tileWidthCRS = Math.abs(
+    forwardTransform(width - 1, 0)[0] - forwardTransform(0, 0)[0],
+  );
+  const roundTripThreshold = tileWidthCRS * 0.01;
+
+  let vi = 0;
+  for (let r = 0; r <= rows; r++) {
+    for (let c = 0; c <= cols; c++) {
+      const u = c / cols;
+      const v = r / rows;
+      const pixelX = u * (width - 1);
+      const pixelY = v * (height - 1);
+      const [ix, iy] = forwardTransform(pixelX, pixelY);
+      const [ox, oy] = forwardReproject(ix, iy);
+      positions[vi * 3] = ox;
+      positions[vi * 3 + 1] = oy;
+      positions[vi * 3 + 2] = 0;
+      texCoords[vi * 2] = u;
+      texCoords[vi * 2 + 1] = v;
+      let ood = !Number.isFinite(ox) || !Number.isFinite(oy);
+      if (!ood && inverseReproject && roundTripThreshold > 0) {
+        const [ix2, iy2] = inverseReproject(ox, oy);
+        if (!Number.isFinite(ix2) || !Number.isFinite(iy2)) {
+          ood = true;
+        } else {
+          const err = Math.sqrt(
+            (ix2 - ix) * (ix2 - ix) + (iy2 - iy) * (iy2 - iy),
+          );
+          if (err > roundTripThreshold) {
+            ood = true;
+          }
+        }
+      }
+      if (ood) {
+        isOOD[vi] = 1;
+      }
+      vi++;
+    }
+  }
+
+  const filteredIndices: number[] = [];
+  for (let r = 0; r < rows; r++) {
+    for (let c = 0; c < cols; c++) {
+      const i0 = r * (cols + 1) + c;
+      const i1 = i0 + 1;
+      const i2 = i0 + (cols + 1);
+      const i3 = i2 + 1;
+      if (!isOOD[i0] && !isOOD[i2] && !isOOD[i1]) {
+        filteredIndices.push(i0, i2, i1);
+      }
+      if (!isOOD[i1] && !isOOD[i2] && !isOOD[i3]) {
+        filteredIndices.push(i1, i2, i3);
+      }
+    }
+  }
+
+  const [positions64Low, positions64High] = splitFloat64Array(positions);
+  const indices = new Uint32Array(filteredIndices);
+  return { indices, positions64High, positions64Low, texCoords };
+}

packages/deck.gl-raster/src/raster-tileset/affine-tileset-level.ts

@@ -44,12 +44,16 @@ export class AffineTilesetLevel implements RasterTilesetLevel {
 
   private readonly _affine: Affine;
   private readonly _invAffine: Affine;
+  private readonly _arrayWidth: number;
+  private readonly _arrayHeight: number;
 
   constructor(options: AffineTilesetLevelOptions) {
     this._affine = options.affine;
     this._invAffine = affine.invert(options.affine);
     this.tileWidth = options.tileWidth;
     this.tileHeight = options.tileHeight;
+    this._arrayWidth = options.arrayWidth;
+    this._arrayHeight = options.arrayHeight;
     this.matrixWidth = Math.ceil(options.arrayWidth / options.tileWidth);
     this.matrixHeight = Math.ceil(options.arrayHeight / options.tileHeight);
 
@@ -78,12 +82,19 @@ export class AffineTilesetLevel implements RasterTilesetLevel {
     const tw = this.tileWidth;
     const th = this.tileHeight;
     const af = this._affine;
-
+    // Clip to actual array extent so corners of the last tile row/column don't
+    // extrapolate past the data boundary. For projections whose valid domain
+    // isn't axis-aligned with the affine (Mollweide, Sinusoidal, Equal Earth),
+    // extrapolated corners fall outside the CRS domain — proj4 maps them to the
+    // pole, collapsing every such corner onto a single ±85.05° Mercator line and
+    // producing a zero-height bounding volume that causes every tile to be culled.
+    const right = Math.min((col + 1) * tw, this._arrayWidth);
+    const bottom = Math.min((row + 1) * th, this._arrayHeight);
     return {
       topLeft: affine.apply(af, col * tw, row * th),
-      topRight: affine.apply(af, (col + 1) * tw, row * th),
-      bottomLeft: affine.apply(af, col * tw, (row + 1) * th),
-      bottomRight: affine.apply(af, (col + 1) * tw, (row + 1) * th),
+      topRight: affine.apply(af, right, row * th),
+      bottomLeft: affine.apply(af, col * tw, bottom),
+      bottomRight: affine.apply(af, right, bottom),
     };
   }
 

packages/deck.gl-raster/src/raster-tileset/raster-tile-traversal.ts

@@ -524,7 +524,7 @@ export class RasterTileNode {
     boundingVolume: OrientedBoundingBox;
     commonSpaceBounds: Bounds;
   } {
-    const [minZ, maxZ] = zRange;
+    const [minZ] = zRange;
 
     const tileCorners = this.level.projectedTileCorners(this.x, this.y);
 
@@ -539,24 +539,38 @@ export class RasterTileNode {
       rescaleEPSG3857ToCommonSpace(xy),
     );
 
-    const refPointPositions: [number, number, number][] = [];
-    for (const p of commonSpacePositions) {
-      refPointPositions.push([p[0], p[1], minZ]);
+    // Filter out NaN positions that arise when tile corners lie outside the
+    // CRS domain (e.g. corners of a bounding rectangle that extend beyond
+    // a Mollweide ellipse). Use only the valid subset for the OBB and AABB.
+    const validPositions = commonSpacePositions.filter(
+      ([x, y]) => Number.isFinite(x) && Number.isFinite(y),
+    );
 
-      if (minZ !== maxZ) {
-        // Also sample at maximum elevation to capture the full 3D volume
-        refPointPositions.push([p[0], p[1], maxZ]);
-      }
+    // If no reference point projects successfully the tile lies entirely
+    // outside the CRS valid domain (e.g. a corner tile whose entire extent
+    // is outside the Mollweide ellipse). Return a degenerate bounding volume
+    // placed far off-screen so both the AABB bounds check and frustum
+    // culling reject it — never selecting it for rendering.
+    if (validPositions.length === 0) {
+      const OFF = -1e8;
+      return {
+        boundingVolume: makeOrientedBoundingBoxFromPoints([
+          [OFF, OFF, minZ],
+          [OFF + 100, OFF, minZ],
+          [OFF, OFF + 100, minZ],
+          [OFF + 100, OFF + 100, minZ],
+        ]),
+        commonSpaceBounds: [OFF, OFF, OFF + 100, OFF + 100],
+      };
     }
 
-    // Compute [minx, miny, maxx, maxy] in common space for quick bounds check
-    // TODO: this doesn't densify edges
+    // Compute [minx, miny, maxx, maxy] in common space for the bounds check
     let minX = Number.POSITIVE_INFINITY;
     let minY = Number.POSITIVE_INFINITY;
     let maxX = Number.NEGATIVE_INFINITY;
     let maxY = Number.NEGATIVE_INFINITY;
 
-    for (const [x, y] of commonSpacePositions) {
+    for (const [x, y] of validPositions) {
       if (x < minX) {
         minX = x;
       }
@@ -572,8 +586,23 @@ export class RasterTileNode {
     }
 
     const commonSpaceBounds: Bounds = [minX, minY, maxX, maxY];
+
+    // Build the OBB from AABB corners with a minimum 1-unit extent in each
+    // axis. makeOrientedBoundingBoxFromPoints produces a degenerate OBB when
+    // the input points are collinear — which happens for projections like
+    // Mollweide where the only valid reference points all fall on the same
+    // latitude (the equatorial band). A degenerate OBB causes
+    // computeVisibility to return -1 (outside frustum), incorrectly culling
+    // the tile. The AABB insideBounds check handles precise culling.
+    const safeMaxX = Math.max(maxX, minX + 1);
+    const safeMaxY = Math.max(maxY, minY + 1);
     return {
-      boundingVolume: makeOrientedBoundingBoxFromPoints(refPointPositions),
+      boundingVolume: makeOrientedBoundingBoxFromPoints([
+        [minX, minY, minZ],
+        [safeMaxX, minY, minZ],
+        [minX, safeMaxY, minZ],
+        [safeMaxX, safeMaxY, minZ],
+      ]),
       commonSpaceBounds,
     };
   }

packages/deck.gl-raster/src/raster-tileset/raster-tileset-2d.ts

@@ -180,8 +180,18 @@ export class RasterTileset2D extends Tileset2D {
     // the tile to keep `RasterLayer`'s reprojection-equality check stable
     // across renders (deck.gl recreates the layer instance every render, so
     // per-render-derived closures would regenerate the mesh every frame).
-    this.projectPosition = (x, y) =>
-      rescaleEPSG3857ToCommonSpace(descriptor.projectTo3857(x, y));
+    this.projectPosition = (x, y) => {
+      const result = rescaleEPSG3857ToCommonSpace(
+        descriptor.projectTo3857(x, y),
+      );
+      if (Number.isFinite(result[0]) && Number.isFinite(result[1])) {
+        return result;
+      }
+      // Point is outside the CRS domain (e.g. beyond the Mollweide ellipse).
+      // Return NaN so the GPU rasterizer discards any triangle that touches
+      // this vertex, producing a clean domain-boundary cutoff with no artifacts.
+      return [NaN, NaN];
+    };
     this.unprojectPosition = (cx, cy) => {
       const [mx, my] = rescaleCommonSpaceToEPSG3857([cx, cy]);
       return descriptor.projectFrom3857(mx, my);

packages/deck.gl-raster/tests/raster-tileset/affine-tileset-level.test.ts

@@ -86,6 +86,37 @@ describe("AffineTilesetLevel", () => {
       expect(corners.bottomRight).toEqual([140, 160]);
     });
 
+    it("clips the right edge of the last-column tile to arrayWidth", () => {
+      // arrayWidth=6 means tile 1 (pixels 4..7) is clipped to pixel 6.
+      const level = new AffineTilesetLevel({
+        affine: SQUARE_AFFINE,
+        arrayWidth: 6,
+        arrayHeight: 8,
+        tileWidth: 4,
+        tileHeight: 4,
+        mpu: 1,
+      });
+      // Interior tile (col=0): right edge is at pixel 4, unclipped.
+      expect(level.projectedTileCorners(0, 0).topRight).toEqual([140, 200]);
+      // Last-column tile (col=1): right edge clipped from pixel 8 to pixel 6.
+      expect(level.projectedTileCorners(1, 0).topRight).toEqual([160, 200]);
+    });
+
+    it("clips the bottom edge of the last-row tile to arrayHeight", () => {
+      const level = new AffineTilesetLevel({
+        affine: SQUARE_AFFINE,
+        arrayWidth: 8,
+        arrayHeight: 6,
+        tileWidth: 4,
+        tileHeight: 4,
+        mpu: 1,
+      });
+      // Interior tile (row=0): bottom at pixel 4, unclipped.
+      expect(level.projectedTileCorners(0, 0).bottomLeft).toEqual([100, 160]);
+      // Last-row tile (row=1): bottom clipped from pixel 8 to pixel 6.
+      expect(level.projectedTileCorners(0, 1).bottomLeft).toEqual([100, 140]);
+    });
+
     it("returns rotated quadrilateral corners for a rotated affine", () => {
       const level = new AffineTilesetLevel({
         affine: ROTATED_AFFINE,