diff --git a/dev-docs/specs/2026-05-27-antimeridian-crossing-tile-design.md b/dev-docs/specs/2026-05-27-antimeridian-crossing-tile-design.md index 8db71e2b..f8e83932 100644 --- a/dev-docs/specs/2026-05-27-antimeridian-crossing-tile-design.md +++ b/dev-docs/specs/2026-05-27-antimeridian-crossing-tile-design.md @@ -33,14 +33,20 @@ projectPosition = (x, y) => rescaleEPSG3857ToCommonSpace(descriptor.projectTo385 Rather than keep the crossing tile as one mesh and fight proj4 to make its coordinates continuous (the **render-as-one** family: #374 output-space shift, #269 reprojector unwrap, #353 global `+over`), **split the tile at the antimeridian into a west piece and an east piece.** Each piece lies wholly on one side of the dateline, so: -- The west piece is monotonic in 3857 (all +x → common-x up to 512); the east piece all −x → common-x from 0. **The discontinuity never exists within a piece.** -- `projectTo3857` stays stock — no unwrap, no proj4 reconfiguration, no `+over`. +- The west piece is monotonic in 3857 (all +x → common-x up to 512); the east piece all −x → common-x from 0. **The discontinuity exists only *at* the shared seam edge, not within a piece's interior** (see "Seam handling" below). +- Almost no projection change: the west piece uses stock `projectTo3857`; the east piece needs only a trivial one-line seam fix — no proj4 reconfiguration, no `+over`, no phase-unwrap. - The `RasterReprojector` needs zero antimeridian awareness — Delatin converges normally on each piece. - Mesh vertices stay within `[0, 512]`, so the fp64 high-zoom precision scheme ([`coordinate-systems.md`](../coordinate-systems.md)) is untouched. - Each piece is a normal tile that the merged world-copy traversal (#518) selects and draws across copies. The antimeridian becomes *a tile boundary*, which the pipeline already handles, instead of a coordinate-space discontinuity. +### Seam handling + +Splitting at the antimeridian is *almost* enough — but not quite. proj4 normalizes ±180° to the **positive** boundary (+max_X / common-x 512). That's correct for the west piece (its right edge *is* +180°), but the east piece's left edge is also the antimeridian and must sit at common-x 0 (−max_X). With stock proj4 the east piece's seam corner lands at 512 while its interior is near 0, so its seed triangle still spans the world and the reprojector diverges — the original #366 failure. + +The fix is local to the **wrapped (negative-side) piece** only: in its `forwardReproject`, **if the projected X comes back positive, subtract one world-width** (the +max boundary → −max). Within that piece the *only* vertex proj4 places on the positive side is the ±180° seam, so this single sign test flips exactly the seam corner and leaves the interior untouched. It is **output-sign-based, not an input-value test** — the seam may be lng +180° or −180° depending on the source's longitude convention, and both pieces share the same seam *input*, so only *which piece* you're rendering decides the handling (known at cut time: the wrapped piece is the one whose interior projects to negative X). `inverseReproject` is unchanged: the piece is now a clean negative range the stock inverse maps back correctly. This is **not** the general phase-unwrap that sank #374 — the piece is known a priori to be wholly on the negative side, so the rule is trivial and deterministic. + ### Why not render-as-one Render-as-one is simpler at the render layer (one mesh, one draw, no internal seam) and more CRS-general (it unwraps the output value, indifferent to source pixel geometry). But it re-attempts the exact unwrap that has failed three times: detection needs phase-unwrapping (a full-world continuous tile must not be mistaken for a crossing tile), mesh vertices leave `[0, 512]` (precision risk), and forward+inverse must stay consistent. We choose cut-in-two as the primary mechanism and keep **render-as-one as the documented fallback for curved-meridian CRS** (see Scope), where cut-in-two degrades. diff --git a/examples/cog-basic/src/App.tsx b/examples/cog-basic/src/App.tsx index 0cee06a9..1db1620b 100644 --- a/examples/cog-basic/src/App.tsx +++ b/examples/cog-basic/src/App.tsx @@ -15,6 +15,13 @@ import type { MapRef } from "react-map-gl/maplibre"; import { Map as MaplibreMap } from "react-map-gl/maplibre"; const COG_OPTIONS: { title: string; url: string; attribution?: ReactNode }[] = [ + { + // Dev-only fixture served by examples/cog-basic/vite.config.ts from the + // geotiff-test-data submodule. EPSG:4326, bbox (−204, −18, −162, 24); + // crosses native −180° at u ≈ 24/42 — exercises antimeridian split. + title: "Antimeridian fixture (dev only)", + url: "/__fixtures/antimeridian.tif", + }, { title: "Sentinel-2 True Color Image (New York, 2026)", url: "https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/18/T/WL/2026/1/S2B_18TWL_20260101_0_L2A/TCI.tif", diff --git a/examples/cog-basic/vite.config.ts b/examples/cog-basic/vite.config.ts index 6b32cedb..8955ea55 100644 --- a/examples/cog-basic/vite.config.ts +++ b/examples/cog-basic/vite.config.ts @@ -1,8 +1,69 @@ +import { createReadStream, promises as fsp } from "node:fs"; +import path from "node:path"; +import { fileURLToPath } from "node:url"; import react from "@vitejs/plugin-react"; +import type { Plugin } from "vite"; import { defineConfig } from "vite"; +const here = path.dirname(fileURLToPath(import.meta.url)); +const fixturesDir = path.resolve( + here, + "../../fixtures/geotiff-test-data/rasterio_generated/fixtures", +); + +/** + * Dev-only middleware that serves files under `examples/cog-basic` at + * `/__fixtures/` from the vendored geotiff-test-data submodule. + * + * Honors HTTP Range requests with 206 + Content-Range so the COGLayer + * behaves identically to a production COG bucket — important for any + * fixture used to validate range-read code paths. + */ +function localFixtures(): Plugin { + return { + name: "local-fixtures", + configureServer(server) { + server.middlewares.use("/__fixtures/", (req, res, next) => { + const requested = decodeURIComponent(req.url ?? "").replace(/^\/+/, ""); + const filePath = path.resolve(fixturesDir, requested); + if ( + filePath !== fixturesDir && + !filePath.startsWith(fixturesDir + path.sep) + ) { + res.statusCode = 403; + res.end(); + return; + } + fsp + .stat(filePath) + .then((stat) => { + res.setHeader("Content-Type", "image/tiff"); + res.setHeader("Accept-Ranges", "bytes"); + const range = req.headers.range; + const m = range ? /^bytes=(\d+)-(\d*)$/.exec(String(range)) : null; + if (m) { + const start = Number.parseInt(m[1]!, 10); + const end = m[2] ? Number.parseInt(m[2], 10) : stat.size - 1; + res.statusCode = 206; + res.setHeader( + "Content-Range", + `bytes ${start}-${end}/${stat.size}`, + ); + res.setHeader("Content-Length", String(end - start + 1)); + createReadStream(filePath, { start, end }).pipe(res); + return; + } + res.setHeader("Content-Length", String(stat.size)); + createReadStream(filePath).pipe(res); + }) + .catch(next); + }); + }, + }; +} + export default defineConfig({ - plugins: [react()], + plugins: [react(), localFixtures()], base: "/deck.gl-raster/examples/cog-basic/", worker: { format: "es" }, server: { diff --git a/packages/deck.gl-raster/src/raster-tile-layer/raster-tile-layer.ts b/packages/deck.gl-raster/src/raster-tile-layer/raster-tile-layer.ts index 8981f843..b6579a3e 100644 --- a/packages/deck.gl-raster/src/raster-tile-layer/raster-tile-layer.ts +++ b/packages/deck.gl-raster/src/raster-tile-layer/raster-tile-layer.ts @@ -1,9 +1,4 @@ -import type { - CompositeLayerProps, - CoordinateSystem, - DefaultProps, - Layer, -} from "@deck.gl/core"; +import type { CompositeLayerProps, DefaultProps, Layer } from "@deck.gl/core"; import { CompositeLayer } from "@deck.gl/core"; import type { _Tile2DHeader as Tile2DHeader, @@ -13,6 +8,7 @@ import type { } from "@deck.gl/geo-layers"; import { TileLayer } from "@deck.gl/geo-layers"; import type { ReprojectionFns } from "@developmentseed/raster-reproject"; +import { triangulateRectangle } from "@developmentseed/raster-reproject"; import type { Device } from "@luma.gl/core"; import { renderDebugTileOutline } from "../layer-utils.js"; import type { RenderTileResult } from "../raster-layer.js"; @@ -380,9 +376,8 @@ export class RasterTileLayer< descriptor: RasterTilesetDescriptor, renderTile: NonNullable["renderTile"]>, ): Layer[] { - const { maxError, debug, debugOpacity } = this.props; + const { debug } = this.props; const tile = props.tile as Tile2DHeader & RasterTileMetadata; - const debugLayers = debug ? this._renderDebug(tile, props.data ?? null) : []; @@ -390,69 +385,155 @@ export class RasterTileLayer< if (!props.data) { return debugLayers; } - - // Access forwardTransform/inverseTransform from tile metadata so that - // reference equality holds across renders. - const { forwardTransform, inverseTransform } = tile; const tileResult = renderTile(props.data); if (!tileResult) { return debugLayers; } - const { image, renderPipeline } = tileResult; - const { width, height } = props.data; const isGlobe = this.context.viewport.resolution !== undefined; - let reprojectionFns: ReprojectionFns; - let coordinateSystem: CoordinateSystem; - if (isGlobe) { - // Globe view - reprojectionFns = { - forwardTransform, - inverseTransform, - forwardReproject: descriptor.projectTo4326, - inverseReproject: descriptor.projectFrom4326, - }; - coordinateSystem = "lnglat"; - } else { - // Web Mercator: render the mesh directly in deck.gl common space. - // - // The tile's `_projectPosition` maps source CRS → common space, support - // high precision with fp64 emulation. - // - // `_projectPosition`/`_unprojectPosition` must be reference-stable across - // renders to avoid regenerating the mesh and recompiling the shader every - // frame. - const { _projectPosition, _unprojectPosition } = tile; - reprojectionFns = { - forwardTransform, - inverseTransform, - forwardReproject: _projectPosition, - inverseReproject: _unprojectPosition, - }; - coordinateSystem = "cartesian"; - } + const rasterLayers = + !isGlobe && tile._antimeridianCut + ? this._renderAntimeridianTile({ + baseId: props.id, + tile, + data: props.data, + tileResult, + uCut: tile._antimeridianCut.uCut, + }) + : this._renderNormalTile({ + baseId: props.id, + tile, + data: props.data, + tileResult, + descriptor, + isGlobe, + }); + return [...rasterLayers, ...debugLayers]; + } - const rasterLayer = new RasterLayer( - this.getSubLayerProps({ - id: `${props.id}-raster`, - width, - height, - // Passing `image: undefined` explicitly would trip isAsyncPropLoading - // and cause a transient black flash (see issue #376). - ...(image !== undefined && { image }), - renderPipeline, - maxError, - reprojectionFns, - // Web Mercator: clamp the mesh to the valid latitude band for tiles - // past ±85.051°. Globe renders the full mesh (it shows the poles). - initialTriangulation: isGlobe - ? undefined - : tile._webMercatorInitialTriangulation, - debug, - debugOpacity, - coordinateSystem, - }), - ); - return [rasterLayer, ...debugLayers]; + /** + * Shared base props for every `RasterLayer` produced from a tile — the + * geometry + render-pipeline pieces that don't depend on globe-vs-mercator + * or whether the tile crosses the antimeridian. + */ + private _baseRasterProps( + data: NonNullable, + tileResult: RenderTileResult, + ) { + const { maxError, debug, debugOpacity } = this.props; + const { image, renderPipeline } = tileResult; + return { + width: data.width, + height: data.height, + // Passing `image: undefined` explicitly would trip isAsyncPropLoading + // and cause a transient black flash (see issue #376). + ...(image !== undefined && { image }), + renderPipeline, + maxError, + debug, + debugOpacity, + }; + } + + /** + * Build the one `RasterLayer` for a tile that renders as a single mesh — + * either the globe path (lng/lat coordinate system, descriptor-level + * projection) or the Web Mercator non-crossing path (cartesian common + * space, per-tile projection with the optional ±85.051° clamp seed). + */ + private _renderNormalTile(opts: { + baseId: string; + tile: Tile2DHeader & RasterTileMetadata; + data: NonNullable; + tileResult: RenderTileResult; + descriptor: RasterTilesetDescriptor; + isGlobe: boolean; + }): Layer[] { + const { baseId, tile, data, tileResult, descriptor, isGlobe } = opts; + const { forwardTransform, inverseTransform } = tile; + + // Web Mercator: render in deck.gl common space using the tile's + // reference-stable `_projectPosition`/`_unprojectPosition` so the mesh + // does not regenerate (and the shader does not recompile) every frame. + // Globe: render in lng/lat using the descriptor's 4326 projection. + const reprojectionFns: ReprojectionFns = isGlobe + ? { + forwardTransform, + inverseTransform, + forwardReproject: descriptor.projectTo4326, + inverseReproject: descriptor.projectFrom4326, + } + : { + forwardTransform, + inverseTransform, + forwardReproject: tile._projectPosition, + inverseReproject: tile._unprojectPosition, + }; + + return [ + new RasterLayer( + this.getSubLayerProps({ + ...this._baseRasterProps(data, tileResult), + id: `${baseId}-raster`, + reprojectionFns, + coordinateSystem: isGlobe ? "lnglat" : "cartesian", + // Globe shows the poles; Web Mercator clamps tiles past ±85.051° + // to the valid latitude band via the seed (undefined when no clamp + // is needed). + initialTriangulation: isGlobe + ? undefined + : tile._webMercatorInitialTriangulation, + }), + ), + ]; + } + + /** + * Build the two `RasterLayer`s for a Web-Mercator tile that crosses ±180°: + * a west piece (UV `[0, uCut]`) and an east piece (UV `[uCut, 1]`). Each + * piece uses its own `ReprojectionFns` bundle from the tile metadata — + * the bundle composes a `+k·360°` longitude shift into the geotransform + * so the piece's native lngs stay inside proj4's valid range, and pairs + * it with the stock `_projectPosition`/`_unprojectPosition` so the + * forward/inverse round-trip cleanly. The two pieces thus render in + * different world copies; deck.gl `repeat: true` + world-copy traversal + * (#518) bring them together visually. The split itself lives in each + * piece's `triangulateRectangle` seed. + */ + private _renderAntimeridianTile(opts: { + baseId: string; + tile: Tile2DHeader & RasterTileMetadata; + data: NonNullable; + tileResult: RenderTileResult; + uCut: number; + }): Layer[] { + const { baseId, tile, data, tileResult, uCut } = opts; + // `antimeridianCut` returns the seam location as a fraction of the tile's + // geographic span (0..1 over the full west→east lng range). `RasterLayer` + // constructs its reprojector with `width + 1` (raster-layer.ts:252), so + // the reprojector's UV*(W-1) = pixel-index maps UV `(seamCol / W)` + // exactly onto the seam pixel — no scaling needed here. + const baseProps = { + ...this._baseRasterProps(data, tileResult), + coordinateSystem: "cartesian" as const, + }; + return [ + new RasterLayer( + this.getSubLayerProps({ + ...baseProps, + id: `${baseId}-raster-west`, + reprojectionFns: tile._westReprojection!, + initialTriangulation: triangulateRectangle(0, 0, uCut, 1), + }), + ), + new RasterLayer( + this.getSubLayerProps({ + ...baseProps, + id: `${baseId}-raster-east`, + reprojectionFns: tile._eastReprojection!, + initialTriangulation: triangulateRectangle(uCut, 0, 1, 1), + }), + ), + ]; } } diff --git a/packages/deck.gl-raster/src/raster-tileset/antimeridian-cut.ts b/packages/deck.gl-raster/src/raster-tileset/antimeridian-cut.ts new file mode 100644 index 00000000..ccd6c348 --- /dev/null +++ b/packages/deck.gl-raster/src/raster-tileset/antimeridian-cut.ts @@ -0,0 +1,85 @@ +/** + * WGS84 longitudes of a tile's four corners, as returned by + * `descriptor.projectTo4326(corner)[0]` — i.e. native and **not** normalized to + * (−180, 180]. For a north-up geotransform the west edge's lng is strictly + * less than the east edge's lng (proj4 4326→4326 is identity and does not + * wrap), so a tile crossing the antimeridian shows up as a span like + * (−204, −162) rather than (156, −162). + */ +export interface CornerLongitudes { + topLeft: number; + topRight: number; + bottomLeft: number; + bottomRight: number; +} + +/** A vertical antimeridian cut in a tile's UV space. */ +export interface AntimeridianCut { + /** UV u-coordinate (0..1) where the tile crosses ±180°. */ + uCut: number; +} + +/** + * Tolerance for treating the top- and bottom-edge crossings as the same u + * (i.e. the cut is vertical). Crossings further apart than this are rejected as + * slanted. + */ +const U_EPSILON = 1e-6; + +/** + * Locate where a single horizontal edge crosses the antimeridian, as a fraction + * of the edge's eastward span (0 at the west corner, 1 at the east corner). + * + * Returns `undefined` if the edge does not cross. Works on native un-normalized + * longitudes (e.g. `westLng = −204`, `eastLng = −162`) by searching for the + * smallest antimeridian line `−180 + 360k` strictly interior to `(westLng, + * eastLng)`. Strict inequalities give the correct non-crossing answer when a + * corner lies exactly on ±180. + */ +function edgeUCut(westLng: number, eastLng: number): number | undefined { + // Degenerate or non-monotonic edge — caller is expected to pass + // west-then-east in the source CRS's native ordering. + if (eastLng <= westLng) { + return undefined; + } + // Smallest antimeridian line (−180 + 360k) strictly greater than westLng. + const k = Math.ceil((westLng + 180) / 360); + const seam = -180 + 360 * k; + if (seam <= westLng || seam >= eastLng) { + return undefined; + } + return (seam - westLng) / (eastLng - westLng); +} + +/** + * Detect whether a tile crosses the antimeridian and, if so, locate the cut. + * + * Only **axis-aligned (vertical) crossings** are handled today (MVP): the top + * and bottom edges must cross ±180° at the same u. We *should* eventually + * support the general case — slanted cuts (rotated geotransforms) and curved + * cuts (non-geographic CRSs) — but for now those return `undefined` and fall + * back to a single full-mesh layer. See issue #575. + * + * Assumes u increases eastward (standard north-up geotransform) and that + * corner longitudes are passed in native, un-normalized form — that is what + * `descriptor.projectTo4326` returns for an EPSG:4326 source whose + * `ModelTiepoint` sits past ±180° (e.g. `−204°`). + */ +export function antimeridianCut( + cornerLngs: CornerLongitudes, +): AntimeridianCut | undefined { + const { topLeft, topRight, bottomLeft, bottomRight } = cornerLngs; + + const topUCut = edgeUCut(topLeft, topRight); + const bottomUCut = edgeUCut(bottomLeft, bottomRight); + if (topUCut === undefined || bottomUCut === undefined) { + return undefined; + } + // Vertical only for now: both edges must cross at the same u. A slanted cut + // (top and bottom crossing at different u) is a valid antimeridian crossing + // we don't yet handle — see the function docstring and issue #575. + if (Math.abs(topUCut - bottomUCut) > U_EPSILON) { + return undefined; + } + return { uCut: (topUCut + bottomUCut) / 2 }; +} diff --git a/packages/deck.gl-raster/src/raster-tileset/raster-tile-traversal.ts b/packages/deck.gl-raster/src/raster-tileset/raster-tile-traversal.ts index fac63439..da286d8e 100644 --- a/packages/deck.gl-raster/src/raster-tileset/raster-tile-traversal.ts +++ b/packages/deck.gl-raster/src/raster-tileset/raster-tile-traversal.ts @@ -54,7 +54,7 @@ import type { * The origin (0,0) is at the top-left corner, and (512,512) is at the * bottom-right. */ -const TILE_SIZE = 512; +export const TILE_SIZE = 512; /** * Maximum number of world copies to test on each side of the primary world diff --git a/packages/deck.gl-raster/src/raster-tileset/raster-tileset-2d.ts b/packages/deck.gl-raster/src/raster-tileset/raster-tileset-2d.ts index f99cb7a8..26378dec 100644 --- a/packages/deck.gl-raster/src/raster-tileset/raster-tileset-2d.ts +++ b/packages/deck.gl-raster/src/raster-tileset/raster-tileset-2d.ts @@ -14,8 +14,13 @@ import type { } from "@deck.gl/geo-layers"; import { _Tileset2D as Tileset2D } from "@deck.gl/geo-layers"; import { transformBounds } from "@developmentseed/proj"; -import type { InitialTriangulation } from "@developmentseed/raster-reproject"; +import type { + InitialTriangulation, + ReprojectionFns, +} from "@developmentseed/raster-reproject"; import type { Matrix4 } from "@math.gl/core"; +import type { AntimeridianCut } from "./antimeridian-cut.js"; +import { antimeridianCut } from "./antimeridian-cut.js"; import { BoundingVolumeCache } from "./bounding-volume-cache.js"; import { getTileIndices, @@ -109,6 +114,32 @@ export type RasterTileMetadata = { * full mesh. See {@link createInitialWebMercatorTriangulation}. */ _webMercatorInitialTriangulation?: InitialTriangulation; + + /** + * Vertical cut at which this tile crosses ±180°, or `undefined` if the tile + * does not cross the antimeridian (or crosses with a slanted/curved cut that + * the MVP does not yet handle). Consumed by `RasterTileLayer._renderSubLayers` + * in the Web Mercator branch to split the tile into a west + east piece. See + * {@link antimeridianCut}. + */ + _antimeridianCut?: AntimeridianCut; + + /** + * Reprojection bundle for the west piece of an antimeridian-crossing tile, + * or `undefined` for non-crossing tiles. The piece's `forwardTransform` + * composes a `+k·360°` longitude shift onto the original geotransform so + * the piece's native longitudes land inside proj4's valid `(−180°, 180°]` + * range — letting the stock `_projectPosition` / `_unprojectPosition` + * round-trip cleanly (which the reprojector's error metric relies on). + * The visual side effect is that the west piece renders in the world-copy + * where its lngs end up after the shift; world-copy traversal places it + * adjacent to the east piece. Built once in `getTileMetadata` for + * reference stability across renders. + */ + _westReprojection?: ReprojectionFns; + + /** East piece counterpart of {@link RasterTileMetadata._westReprojection}. */ + _eastReprojection?: ReprojectionFns; }; /** @@ -396,6 +427,41 @@ export class RasterTileset2D extends Tileset2D { bottomRight: cornerLat(bottomRight), }); + // Detect whether this tile crosses ±180° and locate the vertical cut. + // Corner longitudes are native (as proj4 returns them — un-normalized for a + // 4326 source with an origin past ±180°). See {@link antimeridianCut}. + const cornerLng = (corner: [number, number]) => + this.descriptor.projectTo4326(corner[0], corner[1])[0]; + const tileWestLng = cornerLng(topLeft); + const tileEastLng = cornerLng(topRight); + const _antimeridianCut = antimeridianCut({ + topLeft: tileWestLng, + topRight: tileEastLng, + bottomLeft: cornerLng(bottomLeft), + bottomRight: cornerLng(bottomRight), + }); + + // For each piece of a crossing tile, compose a `+k·360°` longitude shift + // into the geotransform so the piece's native lngs sit inside proj4's + // valid range. The reprojector's error metric uses `inverseReproject` + // round-trip, which only works when proj4 doesn't have to normalize. + let _westReprojection: ReprojectionFns | undefined; + let _eastReprojection: ReprojectionFns | undefined; + if (_antimeridianCut) { + const { uCut } = _antimeridianCut; + const lngAtCut = tileWestLng + uCut * (tileEastLng - tileWestLng); + _westReprojection = this.buildPieceReprojection( + forwardTransform, + inverseTransform, + (tileWestLng + lngAtCut) / 2, + ); + _eastReprojection = this.buildPieceReprojection( + forwardTransform, + inverseTransform, + (lngAtCut + tileEastLng) / 2, + ); + } + return { bbox: { west, @@ -417,6 +483,41 @@ export class RasterTileset2D extends Tileset2D { _projectPosition: this.projectPosition, _unprojectPosition: this.unprojectPosition, _webMercatorInitialTriangulation, + _antimeridianCut, + _westReprojection, + _eastReprojection, + }; + } + + /** + * Build a per-piece reprojection bundle for an antimeridian-crossing tile. + * Picks the `k·360°` longitude shift that brings the piece's native lngs + * (identified by `pieceMidLng`) into proj4's valid range, composes that + * shift into the geotransform, and pairs it with the stock projection + * pair. The composed closures are stable for the tile's lifetime. + */ + private buildPieceReprojection( + forwardTransform: ProjectionFunction, + inverseTransform: ProjectionFunction, + pieceMidLng: number, + ): ReprojectionFns { + const lngShift = -Math.round(pieceMidLng / 360) * 360; + if (lngShift === 0) { + return { + forwardTransform, + inverseTransform, + forwardReproject: this.projectPosition, + inverseReproject: this.unprojectPosition, + }; + } + return { + forwardTransform: (px, py) => { + const [x, y] = forwardTransform(px, py); + return [x + lngShift, y]; + }, + inverseTransform: (x, y) => inverseTransform(x - lngShift, y), + forwardReproject: this.projectPosition, + inverseReproject: this.unprojectPosition, }; } } diff --git a/packages/deck.gl-raster/tests/raster-tileset/antimeridian-cut.test.ts b/packages/deck.gl-raster/tests/raster-tileset/antimeridian-cut.test.ts new file mode 100644 index 00000000..4c040260 --- /dev/null +++ b/packages/deck.gl-raster/tests/raster-tileset/antimeridian-cut.test.ts @@ -0,0 +1,78 @@ +import { describe, expect, it } from "vitest"; +import { antimeridianCut } from "../../src/raster-tileset/antimeridian-cut.js"; + +// cornerLngs are WGS84 longitudes as returned by +// descriptor.projectTo4326(corner)[0] — native and NOT normalized to +// (−180, 180]. For a north-up geotransform, west < east always (proj4 +// 4326→4326 is identity), so a crossing tile shows up as e.g. (−204, −162), +// not (156, −162). +describe("antimeridianCut", () => { + it("returns undefined for a non-crossing tile (west < east, no seam inside)", () => { + expect( + antimeridianCut({ + topLeft: 10, + topRight: 20, + bottomLeft: 10, + bottomRight: 20, + }), + ).toBeUndefined(); + }); + + it("locates a vertical cut for an axis-aligned crossing tile (native un-normalized lngs)", () => { + // antimeridian.tif: ModelTiepoint origin lng −204°, east edge −162°. + const cut = antimeridianCut({ + topLeft: -204, + topRight: -162, + bottomLeft: -204, + bottomRight: -162, + }); + expect(cut).toBeDefined(); + // The −180° seam is at (−180 − (−204)) / (−162 − (−204)) = 24/42 of the + // edge span. + expect(cut?.uCut).toBeCloseTo(24 / 42, 9); + }); + + it("locates the cut equivalently for an in-range crossing span (170, 190)", () => { + const cut = antimeridianCut({ + topLeft: 170, + topRight: 190, + bottomLeft: 170, + bottomRight: 190, + }); + expect(cut?.uCut).toBeCloseTo(0.5, 9); + }); + + it("returns undefined for a slanted (non-vertical) crossing cut", () => { + // Top edge crosses at u = 24/42 ≈ 0.571; bottom edge at u = 10/40 = 0.25. + expect( + antimeridianCut({ + topLeft: -204, + topRight: -162, + bottomLeft: -190, + bottomRight: -150, + }), + ).toBeUndefined(); + }); + + it("returns undefined when a corner lies exactly on the antimeridian (boundary, non-crossing)", () => { + // Strict-inequality seam-finding treats edges that *touch* ±180 as + // non-crossing (the seam is not strictly interior), avoiding a degenerate + // uCut of 0 or 1. + expect( + antimeridianCut({ + topLeft: 160, + topRight: 180, + bottomLeft: 160, + bottomRight: 180, + }), + ).toBeUndefined(); + expect( + antimeridianCut({ + topLeft: -180, + topRight: -160, + bottomLeft: -180, + bottomRight: -160, + }), + ).toBeUndefined(); + }); +}); diff --git a/packages/deck.gl-raster/tests/raster-tileset/raster-tileset-2d-antimeridian.test.ts b/packages/deck.gl-raster/tests/raster-tileset/raster-tileset-2d-antimeridian.test.ts new file mode 100644 index 00000000..e6be8f9f --- /dev/null +++ b/packages/deck.gl-raster/tests/raster-tileset/raster-tileset-2d-antimeridian.test.ts @@ -0,0 +1,149 @@ +import type { _Tileset2DProps as Tileset2DProps } from "@deck.gl/geo-layers"; +import { compose, scale, translation } from "@developmentseed/affine"; +import { describe, expect, it } from "vitest"; +import { AffineTileset } from "../../src/raster-tileset/affine-tileset.js"; +import { AffineTilesetLevel } from "../../src/raster-tileset/affine-tileset-level.js"; +import { RasterTileset2D } from "../../src/raster-tileset/raster-tileset-2d.js"; + +const identity = (x: number, y: number): [number, number] => [x, y]; + +const PROJECTIONS = { + projectTo3857: identity, + projectFrom3857: identity, + projectTo4326: identity, + projectFrom4326: identity, +}; + +function tilesetProps(): Tileset2DProps { + return { getTileData: () => new Promise(() => {}) } as Tileset2DProps; +} + +describe("RasterTileset2D.getTileMetadata — _antimeridianCut", () => { + it("returns _antimeridianCut on a tile whose native lngs cross ±180 (antimeridian.tif shape)", () => { + // antimeridian.tif: rasterio.from_origin(-204, 24, 1, 1), 42×42 EPSG:4326. + // One tile covers the whole image, with native lngs (−204, −162) crossing + // −180° at u = 24/42. + const level = new AffineTilesetLevel({ + affine: compose(translation(-204, 24), scale(1, -1)), + arrayWidth: 42, + arrayHeight: 42, + tileWidth: 42, + tileHeight: 42, + mpu: 1, + }); + const descriptor = new AffineTileset({ levels: [level], ...PROJECTIONS }); + const tileset = new RasterTileset2D(tilesetProps(), descriptor); + + const metadata = tileset.getTileMetadata({ x: 0, y: 0, z: 0 }); + + expect(metadata._antimeridianCut).toBeDefined(); + expect(metadata._antimeridianCut?.uCut).toBeCloseTo(24 / 42, 9); + }); + + it("does NOT set _antimeridianCut on a non-crossing tile", () => { + // A tile entirely east of the antimeridian: native lngs (0, 170). + const level = new AffineTilesetLevel({ + affine: compose(translation(0, 90), scale(1, -1)), + arrayWidth: 170, + arrayHeight: 180, + tileWidth: 170, + tileHeight: 180, + mpu: 1, + }); + const descriptor = new AffineTileset({ levels: [level], ...PROJECTIONS }); + const tileset = new RasterTileset2D(tilesetProps(), descriptor); + + const metadata = tileset.getTileMetadata({ x: 0, y: 0, z: 0 }); + + expect(metadata._antimeridianCut).toBeUndefined(); + expect(metadata._westReprojection).toBeUndefined(); + expect(metadata._eastReprojection).toBeUndefined(); + }); + + it("shifts the west piece's geotransform by +360° so its native lngs round-trip through proj4", () => { + // antimeridian.tif shape again. West piece native lngs (−204°, −180°) are + // outside proj4's valid range. The piece's forwardTransform must add 360° + // to land in [+156°, +180°]; the matching inverseTransform must subtract + // 360° so pixel coords round-trip cleanly. + const level = new AffineTilesetLevel({ + affine: compose(translation(-204, 24), scale(1, -1)), + arrayWidth: 42, + arrayHeight: 42, + tileWidth: 42, + tileHeight: 42, + mpu: 1, + }); + const descriptor = new AffineTileset({ levels: [level], ...PROJECTIONS }); + const tileset = new RasterTileset2D(tilesetProps(), descriptor); + + const metadata = tileset.getTileMetadata({ x: 0, y: 0, z: 0 }); + + expect(metadata._westReprojection).toBeDefined(); + const { forwardTransform: westFwd, inverseTransform: westInv } = + metadata._westReprojection!; + + // Pixel (0, 0) → native (−204, 24) → shifted (+156, 24) + const [wx0, wy0] = westFwd(0, 0); + expect(wx0).toBeCloseTo(156, 9); + expect(wy0).toBeCloseTo(24, 9); + + // Inverse round-trip: shifted (+156, 24) → pixel (0, 0) + const [wpx0, wpy0] = westInv(156, 24); + expect(wpx0).toBeCloseTo(0, 9); + expect(wpy0).toBeCloseTo(0, 9); + }); + + it("leaves the east piece's geotransform unshifted (its lngs already in range)", () => { + const level = new AffineTilesetLevel({ + affine: compose(translation(-204, 24), scale(1, -1)), + arrayWidth: 42, + arrayHeight: 42, + tileWidth: 42, + tileHeight: 42, + mpu: 1, + }); + const descriptor = new AffineTileset({ levels: [level], ...PROJECTIONS }); + const tileset = new RasterTileset2D(tilesetProps(), descriptor); + + const metadata = tileset.getTileMetadata({ x: 0, y: 0, z: 0 }); + + expect(metadata._eastReprojection).toBeDefined(); + const { forwardTransform: eastFwd } = metadata._eastReprojection!; + // East piece native lng midpoint = −171° ∈ [−180°, 180°] → no shift. + // Seam pixel (col 24, row 0) → native (−180, 24), unchanged. + const [ex0, ey0] = eastFwd(24, 0); + expect(ex0).toBeCloseTo(-180, 9); + expect(ey0).toBeCloseTo(24, 9); + }); + + it("shifts the east piece by −360° for a fixture crossing native +180° (e.g. native [170, 190])", () => { + // Mirror of the antimeridian.tif case: tile native lngs (170, 190). East + // piece (180, 190) is outside proj4's range and must shift by −360°. + const level = new AffineTilesetLevel({ + affine: compose(translation(170, 10), scale(1, -1)), + arrayWidth: 20, + arrayHeight: 20, + tileWidth: 20, + tileHeight: 20, + mpu: 1, + }); + const descriptor = new AffineTileset({ levels: [level], ...PROJECTIONS }); + const tileset = new RasterTileset2D(tilesetProps(), descriptor); + + const metadata = tileset.getTileMetadata({ x: 0, y: 0, z: 0 }); + + expect(metadata._antimeridianCut?.uCut).toBeCloseTo(0.5, 9); + expect(metadata._eastReprojection).toBeDefined(); + const { forwardTransform: eastFwd, inverseTransform: eastInv } = + metadata._eastReprojection!; + // East piece native (180, 190) − 360° = (−180, −170). Right pixel col 20 → + // native 190 → shifted −170. + const [ex, ey] = eastFwd(20, 0); + expect(ex).toBeCloseTo(-170, 9); + expect(ey).toBeCloseTo(10, 9); + // Round-trip. + const [epx, epy] = eastInv(-170, 10); + expect(epx).toBeCloseTo(20, 9); + expect(epy).toBeCloseTo(0, 9); + }); +});