feat: Add polar projection support to RasterReprojector

packages/raster-reproject/src/delatin.ts

@@ -70,9 +70,35 @@ export class RasterReprojector {
 
   /**
    * XY Positions in output CRS, computed via exact forward reprojection.
+   *
+   * When the tile crosses the antimeridian, longitude values are normalized
+   * to a continuous range (e.g., [170, 190] instead of [170, -170]) by
+   * applying `_lngOffset` to negative longitudes.
    */
   exactOutputPositions: number[];
 
+  /**
+   * Whether this tile's output positions cross the antimeridian (±180°).
+   * When true, longitudes in `exactOutputPositions` have been shifted
+   * to maintain continuity.
+   */
+  crossesAntimeridian: boolean = false;
+
+  /**
+   * Whether this tile contains or is very near a geographic pole.
+   * When true, the tile's WGS84 output positions have extreme latitude
+   * values (|lat| > 85°) and the longitude range is very wide, indicating
+   * the pole-centered geometry that requires more mesh refinement.
+   */
+  containsPole: boolean = false;
+
+  /**
+   * Longitude offset applied to normalize antimeridian-crossing tiles.
+   * When `crossesAntimeridian` is true, this is 360 (negative longitudes
+   * are shifted by +360). Otherwise 0.
+   */
+  private _lngOffset: number = 0;
+
   /**
    * triangle vertex indices
    */
@@ -127,19 +153,42 @@ export class RasterReprojector {
     const p2 = this._addPoint(0, v1);
     const p3 = this._addPoint(u1, v1);
 
+    // Detect antimeridian crossing: if the longitude range of the 4 corner
+    // vertices exceeds 180°, the tile crosses the ±180° meridian.
+    this._detectAntimeridian();
+
+    // Detect polar tiles: vertices near ±90° latitude with wide longitude
+    // spread indicate a tile at or near a geographic pole.
+    this._detectPole();
+
     // add initial two triangles
     const t0 = this._addTriangle(p3, p0, p2, -1, -1, -1);
     this._addTriangle(p0, p3, p1, t0, -1, -1);
     this._flush();
   }
 
-  // refine the mesh until its maximum error gets below the given one
-  run(maxError: number = DEFAULT_MAX_ERROR): void {
+  /**
+   * Refine the mesh until its maximum error gets below the given one.
+   *
+   * @param maxError - Maximum allowed reprojection error in pixels.
+   * @param opts.maxTriangles - Safety cap on triangle count. Refinement
+   *   stops when this limit is reached even if maxError hasn't been met.
+   *   Useful for polar tiles where convergence is slow due to extreme
+   *   longitude variation near the pole.
+   */
+  run(
+    maxError: number = DEFAULT_MAX_ERROR,
+    opts?: { maxTriangles?: number },
+  ): void {
     if (maxError <= 0) {
       throw new Error("maxError must be positive");
     }
 
-    while (this.getMaxError() > maxError) {
+    const maxTriangles = opts?.maxTriangles ?? Infinity;
+    while (
+      this.getMaxError() > maxError &&
+      this.triangles.length / 3 < maxTriangles
+    ) {
       this.refine();
     }
   }
@@ -164,6 +213,60 @@ export class RasterReprojector {
     this._pendingLen = 0;
   }
 
+  /**
+   * Detect antimeridian crossing from the initial 4 corner vertices.
+   *
+   * If the longitude range exceeds 180°, the tile crosses the antimeridian.
+   * In that case, shift all negative longitudes by +360 to create a
+   * continuous range (e.g., [170, 190] instead of [170, -170]).
+   *
+   * deck.gl handles longitudes outside [-180, 180] correctly.
+   */
+  private _detectAntimeridian(): void {
+    // Check longitude range of the 4 initial vertices
+    let minLng = Infinity;
+    let maxLng = -Infinity;
+    for (let i = 0; i < 4; i++) {
+      const lng = this.exactOutputPositions[i * 2]!;
+      if (lng < minLng) minLng = lng;
+      if (lng > maxLng) maxLng = lng;
+    }
+
+    if (maxLng - minLng > 180) {
+      this.crossesAntimeridian = true;
+      this._lngOffset = 360;
+
+      // Retroactively fix the initial 4 vertices
+      for (let i = 0; i < 4; i++) {
+        const lng = this.exactOutputPositions[i * 2]!;
+        if (lng < 0) {
+          this.exactOutputPositions[i * 2] = lng + 360;
+        }
+      }
+    }
+  }
+
+  /**
+   * Detect whether this tile contains or is very near a geographic pole.
+   *
+   * A tile is considered polar if it crosses the antimeridian (wide
+   * longitude spread) and any vertex has |latitude| > 75°. The
+   * combination of these conditions uniquely identifies tiles in polar
+   * projections — non-polar antimeridian tiles (e.g., UTM zone 1) have
+   * much narrower longitude ranges that don't exceed 180°.
+   */
+  private _detectPole(): void {
+    if (!this.crossesAntimeridian) return;
+
+    for (let i = 0; i < 4; i++) {
+      const lat = this.exactOutputPositions[i * 2 + 1]!;
+      if (Math.abs(lat) > 75) {
+        this.containsPole = true;
+        return;
+      }
+    }
+  }
+
   /**
    * Conversion of upstream's `_findCandidate` for reprojection error handling.
    *
@@ -248,6 +351,10 @@ export class RasterReprojector {
       // Reproject these linearly-interpolated coordinates **from target CRS
       // to input CRS**. This gives us the **exact position in input space**
       // of the linearly interpolated sample point in output space.
+      //
+      // When the tile crosses the antimeridian, the interpolated longitude
+      // may be >180° due to the longitude offset. proj4 handles extended
+      // longitudes correctly, so we pass them through directly.
       const inputCRSSampled = this.reprojectors.inverseReproject(
         outSampleX,
         outSampleY,
@@ -353,10 +460,14 @@ export class RasterReprojector {
       inputPosition[0],
       inputPosition[1],
     );
-    this.exactOutputPositions.push(
-      exactOutputPosition[0]!,
-      exactOutputPosition[1]!,
-    );
+
+    let lng = exactOutputPosition[0]!;
+    // Normalize longitude for antimeridian-crossing tiles
+    if (this._lngOffset !== 0 && lng < 0) {
+      lng += this._lngOffset;
+    }
+
+    this.exactOutputPositions.push(lng, exactOutputPosition[1]!);
 
     return i;
   }

packages/raster-reproject/tests/antimeridian.test.ts

@@ -0,0 +1,120 @@
+import { describe, expect, it } from "vitest";
+import type { ReprojectionFns } from "../src/delatin";
+import { RasterReprojector } from "../src/delatin";
+
+/** Wrap a longitude to [-180, 180] */
+function wrapLng(lng: number): number {
+  return ((((lng + 180) % 360) + 360) % 360) - 180;
+}
+
+/**
+ * Create ReprojectionFns that simulate a raster tile in a source CRS where
+ * the forward transform produces coordinates in source CRS space, and
+ * forwardReproject converts them to WGS84 with longitude wrapping (as proj4
+ * does for geographic output CRS).
+ */
+function makeReprojectionFns(
+  originX: number,
+  originY: number,
+  pixelSizeX: number,
+  pixelSizeY: number,
+  opts?: { wrapLongitude?: boolean },
+): ReprojectionFns {
+  const wrap = opts?.wrapLongitude ?? false;
+  return {
+    forwardTransform(pixelX: number, pixelY: number): [number, number] {
+      return [originX + pixelX * pixelSizeX, originY + pixelY * pixelSizeY];
+    },
+    inverseTransform(crsX: number, crsY: number): [number, number] {
+      return [(crsX - originX) / pixelSizeX, (crsY - originY) / pixelSizeY];
+    },
+    forwardReproject(x: number, y: number): [number, number] {
+      // Simulate proj4 behavior: wrap longitude to [-180, 180]
+      return wrap ? [wrapLng(x), y] : [x, y];
+    },
+    inverseReproject(x: number, y: number): [number, number] {
+      // Identity inverse — extended longitudes (>180) pass through
+      // unchanged, matching how proj4 handles inverse projection
+      // (proj4 normalizes longitude internally before inverse-projecting)
+      return [x, y];
+    },
+  };
+}
+
+describe("antimeridian detection", () => {
+  it("does not flag a tile that does not cross the antimeridian", () => {
+    // A tile centered around longitude 0, latitude 45
+    const fns = makeReprojectionFns(-10, 40, 0.1, -0.05, {
+      wrapLongitude: true,
+    });
+    const reprojector = new RasterReprojector(fns, 200, 200);
+    expect(reprojector.crossesAntimeridian).toBe(false);
+  });
+
+  it("flags a tile that crosses the antimeridian", () => {
+    // Tile spans from 170° to 190° in source CRS.
+    // With wrapLongitude=true, forwardReproject wraps 190° → -170°,
+    // so corner longitudes are [170, -170] — a 340° range that triggers
+    // antimeridian detection.
+    const fns = makeReprojectionFns(170, 40, 0.1, -0.05, {
+      wrapLongitude: true,
+    });
+    const reprojector = new RasterReprojector(fns, 200, 200);
+    expect(reprojector.crossesAntimeridian).toBe(true);
+  });
+
+  it("normalizes longitudes to continuous range when crossing antimeridian", () => {
+    const fns = makeReprojectionFns(170, 40, 0.1, -0.05, {
+      wrapLongitude: true,
+    });
+    const reprojector = new RasterReprojector(fns, 200, 200);
+    reprojector.run(0.5);
+
+    // All longitudes should be in [170, 190] — no jumps to negative values
+    for (let i = 0; i < reprojector.exactOutputPositions.length; i += 2) {
+      const lng = reprojector.exactOutputPositions[i]!;
+      expect(lng).toBeGreaterThanOrEqual(170 - 0.1);
+      expect(lng).toBeLessThanOrEqual(190 + 0.1);
+    }
+  });
+
+  it("mesh triangles do not span more than 180 degrees of longitude", () => {
+    const fns = makeReprojectionFns(170, 40, 0.1, -0.05, {
+      wrapLongitude: true,
+    });
+    const reprojector = new RasterReprojector(fns, 200, 200);
+    reprojector.run(0.5);
+
+    const { triangles, exactOutputPositions } = reprojector;
+    for (let t = 0; t < triangles.length; t += 3) {
+      const a = triangles[t]!;
+      const b = triangles[t + 1]!;
+      const c = triangles[t + 2]!;
+
+      const lngA = exactOutputPositions[a * 2]!;
+      const lngB = exactOutputPositions[b * 2]!;
+      const lngC = exactOutputPositions[c * 2]!;
+
+      const maxLng = Math.max(lngA, lngB, lngC);
+      const minLng = Math.min(lngA, lngB, lngC);
+      expect(maxLng - minLng).toBeLessThan(180);
+    }
+  });
+
+  it("does not affect tiles far from the antimeridian", () => {
+    // A tile centered at longitude 0, no wrapping needed
+    const fns = makeReprojectionFns(-10, 40, 0.1, -0.05, {
+      wrapLongitude: true,
+    });
+    const reprojector = new RasterReprojector(fns, 200, 200);
+    reprojector.run(0.5);
+
+    expect(reprojector.crossesAntimeridian).toBe(false);
+    // All longitudes should be in [-10, 10]
+    for (let i = 0; i < reprojector.exactOutputPositions.length; i += 2) {
+      const lng = reprojector.exactOutputPositions[i]!;
+      expect(lng).toBeGreaterThanOrEqual(-10 - 0.1);
+      expect(lng).toBeLessThanOrEqual(10 + 0.1);
+    }
+  });
+});