From a843354d978cf69fc8e21bc2b014b23fc746daab Mon Sep 17 00:00:00 2001 From: Shane Grigsby Date: Wed, 25 Feb 2026 10:08:52 -0800 Subject: [PATCH 1/4] Add antimeridian crossing support to RasterReprojector MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Detect when a tile's output positions span the ±180° meridian and normalize longitudes to a continuous range (e.g., [170, 190] instead of [170, -170]). This prevents mesh triangles from spanning 340° of longitude and ensures correct GPU-side interpolation for tiles near the date line. --- packages/raster-reproject/src/delatin.ts | 71 +++++++++- .../tests/antimeridian.test.ts | 126 ++++++++++++++++++ 2 files changed, 193 insertions(+), 4 deletions(-) create mode 100644 packages/raster-reproject/tests/antimeridian.test.ts diff --git a/packages/raster-reproject/src/delatin.ts b/packages/raster-reproject/src/delatin.ts index 18588f2a..1df69c65 100644 --- a/packages/raster-reproject/src/delatin.ts +++ b/packages/raster-reproject/src/delatin.ts @@ -70,9 +70,27 @@ 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; + + /** + * 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,6 +145,10 @@ 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(); + // add initial two triangles const t0 = this._addTriangle(p3, p0, p2, -1, -1, -1); this._addTriangle(p0, p3, p1, t0, -1, -1); @@ -164,6 +186,39 @@ 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; + } + } + } + } + /** * Conversion of upstream's `_findCandidate` for reprojection error handling. * @@ -248,6 +303,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 +412,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; } diff --git a/packages/raster-reproject/tests/antimeridian.test.ts b/packages/raster-reproject/tests/antimeridian.test.ts new file mode 100644 index 00000000..a7a8b307 --- /dev/null +++ b/packages/raster-reproject/tests/antimeridian.test.ts @@ -0,0 +1,126 @@ +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); + } + }); +}); From 4adae0be2124ff287ef32772b3f23d3f89693a7b Mon Sep 17 00:00:00 2001 From: Shane Grigsby Date: Thu, 26 Feb 2026 13:15:44 -0800 Subject: [PATCH 2/4] Fix lint issues in antimeridian tests --- packages/raster-reproject/tests/antimeridian.test.ts | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/packages/raster-reproject/tests/antimeridian.test.ts b/packages/raster-reproject/tests/antimeridian.test.ts index a7a8b307..8a69711d 100644 --- a/packages/raster-reproject/tests/antimeridian.test.ts +++ b/packages/raster-reproject/tests/antimeridian.test.ts @@ -23,16 +23,10 @@ function makeReprojectionFns( const wrap = opts?.wrapLongitude ?? false; return { forwardTransform(pixelX: number, pixelY: number): [number, number] { - return [ - originX + pixelX * pixelSizeX, - originY + pixelY * pixelSizeY, - ]; + return [originX + pixelX * pixelSizeX, originY + pixelY * pixelSizeY]; }, inverseTransform(crsX: number, crsY: number): [number, number] { - return [ - (crsX - originX) / pixelSizeX, - (crsY - originY) / pixelSizeY, - ]; + return [(crsX - originX) / pixelSizeX, (crsY - originY) / pixelSizeY]; }, forwardReproject(x: number, y: number): [number, number] { // Simulate proj4 behavior: wrap longitude to [-180, 180] From cac57d168522f5bc5f35ee3647394bdd861bad04 Mon Sep 17 00:00:00 2001 From: Shane Grigsby Date: Wed, 25 Feb 2026 10:17:21 -0800 Subject: [PATCH 3/4] Add polar projection support to RasterReprojector MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Detect tiles containing or near a geographic pole (|lat| > 75° with antimeridian crossing). Add a maxTriangles safety cap to run() to prevent infinite refinement near poles, where the extreme longitude variation can require many more triangles to converge. The existing antimeridian normalization handles the longitude wrapping for polar tiles — the mesh refines around the remaining discontinuity. deck.gl's GlobeView renders the resulting mesh correctly since projectPosition handles arbitrary (lng, lat) values. --- packages/raster-reproject/src/delatin.ts | 54 +++++- packages/raster-reproject/tests/polar.test.ts | 168 ++++++++++++++++++ 2 files changed, 219 insertions(+), 3 deletions(-) create mode 100644 packages/raster-reproject/tests/polar.test.ts diff --git a/packages/raster-reproject/src/delatin.ts b/packages/raster-reproject/src/delatin.ts index 1df69c65..b8b3ec1a 100644 --- a/packages/raster-reproject/src/delatin.ts +++ b/packages/raster-reproject/src/delatin.ts @@ -84,6 +84,14 @@ export class RasterReprojector { */ 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 @@ -149,19 +157,38 @@ export class RasterReprojector { // 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(); } } @@ -219,6 +246,27 @@ export class RasterReprojector { } } + /** + * 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. * diff --git a/packages/raster-reproject/tests/polar.test.ts b/packages/raster-reproject/tests/polar.test.ts new file mode 100644 index 00000000..226d42d7 --- /dev/null +++ b/packages/raster-reproject/tests/polar.test.ts @@ -0,0 +1,168 @@ +import { describe, expect, it } from "vitest"; +import type { ReprojectionFns } from "../src/delatin"; +import { RasterReprojector } from "../src/delatin"; + +const R = 6378137; // WGS84 semi-major axis (meters) + +/** + * Simplified south polar stereographic projection for testing. + * + * Maps a rectangular tile in a planar CRS centered on the South Pole + * to WGS84 (lon, lat). The projection is a true polar stereographic + * with the origin at the South Pole (0, 0) in projected coordinates. + * + * Convention: + * - x axis → 90°E direction + * - y axis → 0° (prime meridian) direction + * - South Pole is at (0, 0) in projected coords = (*, -90°) in WGS84 + */ +function makePolarReprojectionFns( + originX: number, + originY: number, + pixelSizeX: number, + pixelSizeY: number, +): ReprojectionFns { + 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] { + // Polar stereographic → WGS84 + const rho = Math.sqrt(x * x + y * y); + if (rho < 0.01) { + // At the pole — longitude is undefined, return 0 by convention + return [0, -90]; + } + const c = 2 * Math.atan2(rho, 2 * R); + const lat = ((c - Math.PI / 2) * 180) / Math.PI; + const lon = (Math.atan2(x, y) * 180) / Math.PI; + return [lon, lat]; + }, + inverseReproject(lon: number, lat: number): [number, number] { + // WGS84 → polar stereographic + const latRad = (lat * Math.PI) / 180; + const lonRad = (lon * Math.PI) / 180; + const rho = 2 * R * Math.tan(Math.PI / 4 + latRad / 2); + const x = rho * Math.sin(lonRad); + const y = rho * Math.cos(lonRad); + return [x, y]; + }, + }; +} + +describe("polar projection support", () => { + it("detects a tile containing the south pole", () => { + // Tile centered on the south pole: spans ±500km around the pole + // Origin is top-left corner in GeoTIFF convention + const fns = makePolarReprojectionFns(-500000, 500000, 5000, -5000); + const reprojector = new RasterReprojector(fns, 200, 200); + + expect(reprojector.crossesAntimeridian).toBe(true); + expect(reprojector.containsPole).toBe(true); + }); + + it("does not flag a tile far from the pole", () => { + // Tile at ~1500-2000km from pole (roughly -75° to -70° latitude) + // At this distance, longitude spread is limited + const fns = makePolarReprojectionFns(1500000, 2500000, 5000, -5000); + const reprojector = new RasterReprojector(fns, 200, 200); + + expect(reprojector.containsPole).toBe(false); + }); + + it("generates a mesh for a polar tile with maxTriangles cap", () => { + // Tile containing the south pole + const fns = makePolarReprojectionFns(-500000, 500000, 5000, -5000); + const reprojector = new RasterReprojector(fns, 200, 200); + + // Run with a safety cap — polar tiles may need many triangles + reprojector.run(2.0, { maxTriangles: 5000 }); + + const numTriangles = reprojector.triangles.length / 3; + expect(numTriangles).toBeGreaterThan(2); + // Each refine() step adds 2+ triangles atomically, so the final count + // may slightly exceed the cap + expect(numTriangles).toBeLessThanOrEqual(5010); + + // Verify mesh has vertices + expect(reprojector.uvs.length).toBeGreaterThan(8); // more than initial 4 vertices + }); + + it("generates a valid mesh for a near-pole tile", () => { + // Tile near the pole but not containing it + // Offset so the pole is outside the tile extent + const fns = makePolarReprojectionFns(200000, 700000, 5000, -5000); + const reprojector = new RasterReprojector(fns, 200, 200); + + // Should still converge with a reasonable triangle count + reprojector.run(2.0, { maxTriangles: 5000 }); + + const numTriangles = reprojector.triangles.length / 3; + expect(numTriangles).toBeGreaterThan(2); + }); + + it("mesh covers the expected latitude range", () => { + // Tile containing the pole, spanning ~±4.5° from the pole + const fns = makePolarReprojectionFns(-500000, 500000, 5000, -5000); + const reprojector = new RasterReprojector(fns, 200, 200); + reprojector.run(2.0, { maxTriangles: 5000 }); + + // Check that output positions include latitudes near -90° + let minLat = Infinity; + let maxLat = -Infinity; + for (let i = 0; i < reprojector.exactOutputPositions.length; i += 2) { + const lat = reprojector.exactOutputPositions[i + 1]!; + if (lat < minLat) minLat = lat; + if (lat > maxLat) maxLat = lat; + } + + // The tile should reach close to the pole (with capped triangle count, + // the mesh may not have a vertex at exactly -90°) + expect(minLat).toBeLessThan(-86); + // And extend away from the pole + expect(maxLat).toBeGreaterThan(-86); + }); + + it("maxTriangles stops refinement early", () => { + const fns = makePolarReprojectionFns(-500000, 500000, 5000, -5000); + const reprojector = new RasterReprojector(fns, 200, 200); + + // Use a very tight error threshold and low maxTriangles + reprojector.run(0.01, { maxTriangles: 100 }); + + const numTriangles = reprojector.triangles.length / 3; + // Each refine() step adds 2+ triangles atomically, so the final count + // may slightly exceed the cap + expect(numTriangles).toBeLessThanOrEqual(110); + // Error should still be above threshold since we stopped early + expect(reprojector.getMaxError()).toBeGreaterThan(0.01); + }); + + it("round-trips through forward and inverse reprojection", () => { + // Verify our test projection is self-consistent + const fns = makePolarReprojectionFns(-500000, 500000, 5000, -5000); + + // Test a few points + const testPoints: [number, number][] = [ + [100000, 200000], + [-300000, 100000], + [0, -400000], + ]; + + for (const [x, y] of testPoints) { + const [lon, lat] = fns.forwardReproject(x, y); + const [x2, y2] = fns.inverseReproject(lon, lat); + expect(x2).toBeCloseTo(x, 0); + expect(y2).toBeCloseTo(y, 0); + } + }); +}); From ddb38de41960ed9a69dbaa110fdb86e17a8d8737 Mon Sep 17 00:00:00 2001 From: Shane Grigsby Date: Thu, 26 Feb 2026 13:16:19 -0800 Subject: [PATCH 4/4] Fix lint issues in polar tests --- packages/raster-reproject/tests/polar.test.ts | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/packages/raster-reproject/tests/polar.test.ts b/packages/raster-reproject/tests/polar.test.ts index 226d42d7..70cbf702 100644 --- a/packages/raster-reproject/tests/polar.test.ts +++ b/packages/raster-reproject/tests/polar.test.ts @@ -24,16 +24,10 @@ function makePolarReprojectionFns( ): ReprojectionFns { return { forwardTransform(pixelX: number, pixelY: number): [number, number] { - return [ - originX + pixelX * pixelSizeX, - originY + pixelY * pixelSizeY, - ]; + return [originX + pixelX * pixelSizeX, originY + pixelY * pixelSizeY]; }, inverseTransform(crsX: number, crsY: number): [number, number] { - return [ - (crsX - originX) / pixelSizeX, - (crsY - originY) / pixelSizeY, - ]; + return [(crsX - originX) / pixelSizeX, (crsY - originY) / pixelSizeY]; }, forwardReproject(x: number, y: number): [number, number] { // Polar stereographic → WGS84