diff --git a/packages/raster-reproject/src/delatin.ts b/packages/raster-reproject/src/delatin.ts index 18588f2a..b8b3ec1a 100644 --- a/packages/raster-reproject/src/delatin.ts +++ b/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; } diff --git a/packages/raster-reproject/tests/antimeridian.test.ts b/packages/raster-reproject/tests/antimeridian.test.ts new file mode 100644 index 00000000..8a69711d --- /dev/null +++ b/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); + } + }); +}); diff --git a/packages/raster-reproject/tests/polar.test.ts b/packages/raster-reproject/tests/polar.test.ts new file mode 100644 index 00000000..70cbf702 --- /dev/null +++ b/packages/raster-reproject/tests/polar.test.ts @@ -0,0 +1,162 @@ +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); + } + }); +});