From 7280242df1d15176693c9e9f6b21daf9bc7aad30 Mon Sep 17 00:00:00 2001 From: jameswillis Date: Wed, 25 Mar 2026 14:28:39 -0700 Subject: [PATCH] fix: handle antimeridian-crossing tiles in Web Mercator reprojection MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Tiles near ±180° longitude cause RasterReprojector mesh refinement to diverge because forwardTo3857 maps nearby source coordinates on opposite sides of the antimeridian to EPSG:3857 x-values ~40M meters apart. Detect these tiles by checking if corner x-values in EPSG:3857 span more than half the world circumference, then wrap the projection functions to make the coordinate space continuous. Extract the wrapping logic into wrapAntimeridianProjections() in proj.ts with unit tests covering no-op passthrough, wrapping detection, coordinate continuity, and forward/inverse round-tripping. Closes #366 --- packages/deck.gl-geotiff/src/cog-layer.ts | 28 +++++++- packages/deck.gl-geotiff/src/proj.ts | 36 +++++++++++ packages/deck.gl-geotiff/tests/proj.test.ts | 71 ++++++++++++++++++++- 3 files changed, 131 insertions(+), 4 deletions(-) diff --git a/packages/deck.gl-geotiff/src/cog-layer.ts b/packages/deck.gl-geotiff/src/cog-layer.ts index 283c2a68..c9bfb69e 100644 --- a/packages/deck.gl-geotiff/src/cog-layer.ts +++ b/packages/deck.gl-geotiff/src/cog-layer.ts @@ -39,7 +39,11 @@ import type { TextureDataT } from "./geotiff/render-pipeline.js"; import { inferRenderPipeline } from "./geotiff/render-pipeline.js"; import { fromAffine } from "./geotiff-reprojection.js"; import type { EpsgResolver } from "./proj.js"; -import { epsgResolver, makeClampedForwardTo3857 } from "./proj.js"; +import { + epsgResolver, + makeClampedForwardTo3857, + wrapAntimeridianProjections, +} from "./proj.js"; /** Size of deck.gl's common coordinate space in world units. * @@ -463,11 +467,29 @@ export class COGLayer< }; deckProjectionProps = {}; } else { + // Wrap projection fns for antimeridian-crossing tiles so the 3857 + // x-coordinates are continuous and RasterReprojector can converge (#366). + const corners = [ + [0, 0], + [width, 0], + [width, height], + [0, height], + ] as const; + const cornerXs = corners.map(([cx, cy]) => { + const [sx, sy] = forwardTransform(cx, cy); + return forwardTo3857(sx, sy)[0]; + }); + const wrapped = wrapAntimeridianProjections( + cornerXs, + forwardTo3857, + inverseFrom3857, + ); + reprojectionFns = { forwardTransform, inverseTransform, - forwardReproject: forwardTo3857, - inverseReproject: inverseFrom3857, + forwardReproject: wrapped.forwardTo3857, + inverseReproject: wrapped.inverseFrom3857, }; // Scale 3857 meters → deck.gl world units (512×512). // diff --git a/packages/deck.gl-geotiff/src/proj.ts b/packages/deck.gl-geotiff/src/proj.ts index 845c9a71..acd81adb 100644 --- a/packages/deck.gl-geotiff/src/proj.ts +++ b/packages/deck.gl-geotiff/src/proj.ts @@ -8,6 +8,11 @@ const WGS84_ELLIPSOID_A = 6378137; // Beyond this, the Mercator projection is undefined. const MAX_WEB_MERCATOR_LAT = 85.05112877980659; +const WEB_MERCATOR_METER_CIRCUMFERENCE = 2 * Math.PI * WGS84_ELLIPSOID_A; +const HALF_CIRCUMFERENCE = WEB_MERCATOR_METER_CIRCUMFERENCE / 2; + +type ProjectionFn = (x: number, y: number) => [number, number]; + /** * Convert a WGS84 longitude/latitude to EPSG:3857 meters analytically. * Valid for latitudes in [-MAX_WEB_MERCATOR_LAT, MAX_WEB_MERCATOR_LAT]. @@ -19,6 +24,37 @@ function wgs84To3857(lon: number, lat: number): [number, number] { return [x, y]; } +/** + * If a tile's EPSG:3857 corner x-values span more than half the globe, wrap + * `forwardTo3857` / `inverseFrom3857` so the coordinate space is continuous. + * + * Returns the original functions unchanged when no wrapping is needed. + */ +export function wrapAntimeridianProjections( + cornerXs3857: number[], + forwardTo3857: ProjectionFn, + inverseFrom3857: ProjectionFn, +): { forwardTo3857: ProjectionFn; inverseFrom3857: ProjectionFn } { + const xMin = Math.min(...cornerXs3857); + const xMax = Math.max(...cornerXs3857); + + if (xMax - xMin <= HALF_CIRCUMFERENCE) { + return { forwardTo3857, inverseFrom3857 }; + } + + return { + forwardTo3857: (x: number, y: number): [number, number] => { + const [px, py] = forwardTo3857(x, y); + return [px < 0 ? px + WEB_MERCATOR_METER_CIRCUMFERENCE : px, py]; + }, + inverseFrom3857: (x: number, y: number): [number, number] => { + const unwrapped = + x > HALF_CIRCUMFERENCE ? x - WEB_MERCATOR_METER_CIRCUMFERENCE : x; + return inverseFrom3857(unwrapped, y); + }, + }; +} + /** * Wrap a proj4 forward projection to EPSG:3857 so that it never returns NaN. * diff --git a/packages/deck.gl-geotiff/tests/proj.test.ts b/packages/deck.gl-geotiff/tests/proj.test.ts index fe76a06b..3bcefb2b 100644 --- a/packages/deck.gl-geotiff/tests/proj.test.ts +++ b/packages/deck.gl-geotiff/tests/proj.test.ts @@ -1,6 +1,9 @@ import proj4 from "proj4"; import { describe, expect, it } from "vitest"; -import { makeClampedForwardTo3857 } from "../src/proj.js"; +import { + makeClampedForwardTo3857, + wrapAntimeridianProjections, +} from "../src/proj.js"; const WGS84_ELLIPSOID_A = 6378137; const EPSG_3857_HALF_CIRCUMFERENCE = Math.PI * WGS84_ELLIPSOID_A; @@ -44,3 +47,69 @@ describe("makeClampedForwardTo3857", () => { expect(y).toBeCloseTo(EPSG_3857_HALF_CIRCUMFERENCE, 0); }); }); + +describe("wrapAntimeridianProjections", () => { + const converter3857 = proj4("EPSG:4326", "EPSG:3857"); + + const forwardTo3857 = (x: number, y: number): [number, number] => + converter3857.forward([x, y], false); + const inverseFrom3857 = (x: number, y: number): [number, number] => + converter3857.inverse([x, y], false); + + it("returns original functions when tile does not cross antimeridian", () => { + // Tile from lon 10° to 20° — well within one hemisphere + const cornerXs = [10, 20].map((lon) => forwardTo3857(lon, 0)[0]); + const result = wrapAntimeridianProjections( + cornerXs, + forwardTo3857, + inverseFrom3857, + ); + expect(result.forwardTo3857).toBe(forwardTo3857); + expect(result.inverseFrom3857).toBe(inverseFrom3857); + }); + + it("wraps functions when tile crosses the antimeridian", () => { + // Tile corners at lon +170° and -170° (crosses ±180°) + const cornerXs = [170, -170].map((lon) => forwardTo3857(lon, 0)[0]); + const result = wrapAntimeridianProjections( + cornerXs, + forwardTo3857, + inverseFrom3857, + ); + // Should return new (wrapped) functions + expect(result.forwardTo3857).not.toBe(forwardTo3857); + expect(result.inverseFrom3857).not.toBe(inverseFrom3857); + }); + + it("produces continuous x-values for antimeridian-crossing tiles", () => { + const cornerXs = [170, -170].map((lon) => forwardTo3857(lon, 0)[0]); + const { forwardTo3857: wrapped } = wrapAntimeridianProjections( + cornerXs, + forwardTo3857, + inverseFrom3857, + ); + + const x170 = wrapped(170, 0)[0]; + const xNeg170 = wrapped(-170, 0)[0]; + + // Both should now be positive and close together (~20° apart in meters) + expect(x170).toBeGreaterThan(0); + expect(xNeg170).toBeGreaterThan(0); + expect(Math.abs(xNeg170 - x170)).toBeLessThan(5_000_000); + }); + + it("round-trips through wrapped forward and inverse", () => { + const cornerXs = [170, -170].map((lon) => forwardTo3857(lon, 0)[0]); + const result = wrapAntimeridianProjections( + cornerXs, + forwardTo3857, + inverseFrom3857, + ); + + // Forward then inverse should recover the original lon/lat + const [mx, my] = result.forwardTo3857(-175, 45); + const [lon, lat] = result.inverseFrom3857(mx, my); + expect(lon).toBeCloseTo(-175, 4); + expect(lat).toBeCloseTo(45, 4); + }); +});