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); + }); +});