diff --git a/examples/globe-view/index.html b/examples/globe-view/index.html new file mode 100644 index 00000000..3a402509 --- /dev/null +++ b/examples/globe-view/index.html @@ -0,0 +1,22 @@ + + + + + + COGLayer Globe View Example + + + +
+ + + diff --git a/examples/globe-view/package.json b/examples/globe-view/package.json new file mode 100644 index 00000000..3d75e7e8 --- /dev/null +++ b/examples/globe-view/package.json @@ -0,0 +1,32 @@ +{ + "name": "deck.gl-cog-globe-example", + "private": true, + "type": "module", + "scripts": { + "dev": "vite", + "build": "vite build", + "preview": "vite preview" + }, + "dependencies": { + "@deck.gl/core": "^9.2.7", + "@deck.gl/geo-layers": "^9.2.7", + "@deck.gl/layers": "^9.2.7", + "@deck.gl/mesh-layers": "^9.2.7", + "@deck.gl/react": "^9.2.7", + "@developmentseed/geotiff": "workspace:^", + "@developmentseed/deck.gl-geotiff": "workspace:^", + "@developmentseed/deck.gl-raster": "workspace:^", + "@luma.gl/core": "9.2.6", + "@luma.gl/shadertools": "9.2.6", + "@luma.gl/webgl": "9.2.6", + "proj4": "^2.20.2", + "react": "^19.2.4", + "react-dom": "^19.2.4" + }, + "devDependencies": { + "@types/react": "^19.2.10", + "@types/react-dom": "^19.2.3", + "@vitejs/plugin-react": "^5.1.2", + "vite": "^7.3.1" + } +} diff --git a/examples/globe-view/src/App.tsx b/examples/globe-view/src/App.tsx new file mode 100644 index 00000000..991b4f87 --- /dev/null +++ b/examples/globe-view/src/App.tsx @@ -0,0 +1,181 @@ +import { _GlobeView as GlobeView } from "@deck.gl/core"; +import { GeoJsonLayer, SolidPolygonLayer } from "@deck.gl/layers"; +import { DeckGL } from "@deck.gl/react"; +import { COGLayer } from "@developmentseed/deck.gl-geotiff"; +import { luma } from "@luma.gl/core"; +import { webgl2Adapter } from "@luma.gl/webgl"; +import { useCallback, useState } from "react"; + +// Register WebGL adapter — required when DeckGL creates its own context +// (unlike MapboxOverlay which reuses MaplibreGL's existing context) +luma.registerAdapters([webgl2Adapter]); + +// New Zealand imagery (NZTM2000 projection) +const COG_URL = + "https://nz-imagery.s3-ap-southeast-2.amazonaws.com/new-zealand/new-zealand_2024-2025_10m/rgb/2193/CC11.tiff"; + +// Antarctic sea ice (polar stereographic) +// const COG_URL = +// "https://data.source.coop/ausantarctic/ghrsst-mur-v2/2020/12/12/20201212090000-JPL-L4_GHRSST-SSTfnd-MUR-GLOB-v02.0-fv04.1_sea_ice_fraction.tif"; + +export default function App() { + const [debug, setDebug] = useState(false); + const [debugOpacity, setDebugOpacity] = useState(0.25); + + // Use initialViewState (uncontrolled) so deck.gl manages view state internally. + // Updating the object reference triggers deck.gl to transition to the new view. + const [initialViewState, setInitialViewState] = useState({ + longitude: 0, + latitude: 0, + zoom: 1, + }); + + const onGeoTIFFLoad = useCallback( + ( + _tiff: unknown, + options: { + geographicBounds: { + west: number; + south: number; + east: number; + north: number; + }; + }, + ) => { + const { west, south, east, north } = options.geographicBounds; + const lonSpan = east - west; + const latSpan = north - south; + const maxSpan = Math.max(lonSpan, latSpan); + // At zoom N, ~360/2^N degrees are visible; subtract 1 for padding + const zoom = Math.log2(360 / maxSpan) - 1; + setInitialViewState({ + longitude: (west + east) / 2, + latitude: (south + north) / 2, + zoom, + }); + }, + [], + ); + + const layers = [ + // Dark background sphere + new SolidPolygonLayer({ + id: "background", + data: [ + [ + [-180, 90], + [0, 90], + [180, 90], + [180, -90], + [0, -90], + [-180, -90], + ], + ], + getPolygon: (d) => d, + stroked: false, + filled: true, + getFillColor: [10, 20, 40], + }), + // Land masses basemap (Natural Earth via deck.gl CDN) + new GeoJsonLayer({ + id: "basemap", + data: "https://d2ad6b4ur7yvpq.cloudfront.net/naturalearth-3.3.0/ne_50m_admin_0_scale_rank.geojson", + stroked: true, + filled: true, + lineWidthMinPixels: 1, + getLineColor: [40, 60, 90], + getFillColor: [25, 40, 70], + }), + new COGLayer({ + id: "cog-layer", + geotiff: COG_URL, + debug, + debugOpacity, + onGeoTIFFLoad, + }), + ]; + + return ( +
+ + + {/* UI Controls */} +
+

+ COGLayer Globe View +

+

+ Displaying COG imagery on a 3D globe +

+ +
+ + + {debug && ( +
+ +
+ )} +
+
+
+ ); +} diff --git a/examples/globe-view/src/main.tsx b/examples/globe-view/src/main.tsx new file mode 100644 index 00000000..f8fc6f51 --- /dev/null +++ b/examples/globe-view/src/main.tsx @@ -0,0 +1,9 @@ +import { StrictMode } from "react"; +import { createRoot } from "react-dom/client"; +import App from "./App"; + +createRoot(document.getElementById("root")!).render( + + + , +); diff --git a/examples/globe-view/tsconfig.json b/examples/globe-view/tsconfig.json new file mode 100644 index 00000000..f0a23505 --- /dev/null +++ b/examples/globe-view/tsconfig.json @@ -0,0 +1,24 @@ +{ + "compilerOptions": { + "target": "ES2020", + "useDefineForClassFields": true, + "lib": ["ES2020", "DOM", "DOM.Iterable"], + "module": "ESNext", + "skipLibCheck": true, + + /* Bundler mode */ + "moduleResolution": "bundler", + "allowImportingTsExtensions": true, + "isolatedModules": true, + "moduleDetection": "force", + "noEmit": true, + "jsx": "react-jsx", + + /* Linting */ + "strict": true, + "noUnusedLocals": true, + "noUnusedParameters": true, + "noFallthroughCasesInSwitch": true + }, + "include": ["src"] +} diff --git a/examples/globe-view/vite.config.ts b/examples/globe-view/vite.config.ts new file mode 100644 index 00000000..26bf8080 --- /dev/null +++ b/examples/globe-view/vite.config.ts @@ -0,0 +1,10 @@ +import react from "@vitejs/plugin-react"; +import { defineConfig } from "vite"; + +export default defineConfig({ + plugins: [react()], + base: "/deck.gl-raster/examples/globe-view/", + server: { + port: 3001, + }, +}); 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 fbcbe026..78b639dd 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 @@ -16,7 +16,7 @@ */ import type { Viewport } from "@deck.gl/core"; -import { _GlobeViewport, assert } from "@deck.gl/core"; +import { _GlobeViewport } from "@deck.gl/core"; import type { TileMatrix, TileMatrixSet } from "@developmentseed/morecantile"; import { xy_bounds } from "@developmentseed/morecantile"; import type { OrientedBoundingBox } from "@math.gl/culling"; @@ -141,6 +141,7 @@ export class RasterTileNode { private _children?: RasterTileNode[] | null; private projectTo3857: ProjectionFunction; + private projectTo4326: ProjectionFunction; constructor( x: number, @@ -149,13 +150,19 @@ export class RasterTileNode { { metadata, projectTo3857, - }: { metadata: TileMatrixSet; projectTo3857: ProjectionFunction }, + projectTo4326, + }: { + metadata: TileMatrixSet; + projectTo3857: ProjectionFunction; + projectTo4326: ProjectionFunction; + }, ) { this.x = x; this.y = y; this.z = z; this.metadata = metadata; this.projectTo3857 = projectTo3857; + this.projectTo4326 = projectTo4326; } /** Get overview info for this tile's z level */ @@ -198,13 +205,14 @@ export class RasterTileNode { const children: RasterTileNode[] = []; - const { metadata, projectTo3857 } = this; + const { metadata, projectTo3857, projectTo4326 } = this; for (let y = minRow; y <= maxRow; y++) { for (let x = minCol; x <= maxCol; x++) { children.push( new RasterTileNode(x, y, childZ, { metadata, projectTo3857, + projectTo4326, }), ); } @@ -263,10 +271,8 @@ export class RasterTileNode { } = params; // Get bounding volume for this tile - const { boundingVolume, commonSpaceBounds } = this.getBoundingVolume( - elevationBounds, - project, - ); + const { boundingVolume, commonSpaceBounds, centerLatitude } = + this.getBoundingVolume(elevationBounds, project); // Step 1: Bounds checking // If geographic bounds are specified, reject tiles outside those bounds @@ -288,8 +294,8 @@ export class RasterTileNode { // Only select this tile if no child is visible (prevents overlapping tiles) // “When pitch is low, force selection at maxZ.” if (!this.childVisible && this.z >= minZ) { - const metersPerScreenPixel = getMetersPerPixelAtBoundingVolume( - boundingVolume, + const metersPerScreenPixel = getMetersPerPixel( + centerLatitude, viewport.zoom, ); // console.log("metersPerScreenPixel", metersPerScreenPixel); @@ -396,13 +402,15 @@ export class RasterTileNode { getBoundingVolume( zRange: ZRange, project: ((xyz: number[]) => number[]) | null, - ): { boundingVolume: OrientedBoundingBox; commonSpaceBounds: Bounds } { - // Case 1: Globe view - need to construct an oriented bounding box from - // reprojected sample points, but also using the `project` param + ): { + boundingVolume: OrientedBoundingBox; + commonSpaceBounds: Bounds; + centerLatitude: number; + } { + // Case 1: Globe view - construct an oriented bounding box from sample + // points projected into globe common space via viewport.projectPosition if (project) { - assert(false, "TODO: implement getBoundingVolume in Globe view"); - // Reproject positions to wgs84 instead, then pass them into `project` - // return makeOrientedBoundingBoxFromPoints(refPointPositions); + return this._getGlobeBoundingVolume(zRange, project); } // (Future) Case 2: Web Mercator input image, can directly compute AABB in @@ -425,6 +433,7 @@ export class RasterTileNode { private _getGenericBoundingVolume(zRange: ZRange): { boundingVolume: OrientedBoundingBox; commonSpaceBounds: Bounds; + centerLatitude: number; } { const tileMatrix = this.tileMatrix; const [minZ, maxZ] = zRange; @@ -469,9 +478,81 @@ export class RasterTileNode { } const commonSpaceBounds: Bounds = [minX, minY, maxX, maxY]; + const boundingVolume = makeOrientedBoundingBoxFromPoints(refPointPositions); + const [, centerLatitude] = worldToLngLat(boundingVolume.center); + + return { + boundingVolume, + commonSpaceBounds, + centerLatitude, + }; + } + + /** + * Globe view bounding volume. + * + * Sample reference points, reproject to WGS84, then project into globe + * common space via viewport.projectPosition. + */ + private _getGlobeBoundingVolume( + zRange: ZRange, + project: (xyz: number[]) => number[], + ): { + boundingVolume: OrientedBoundingBox; + commonSpaceBounds: Bounds; + centerLatitude: number; + } { + const tileMatrix = this.tileMatrix; + const [minZ, maxZ] = zRange; + + const tileCrsBounds = computeProjectedTileBounds(tileMatrix, { + x: this.x, + y: this.y, + }); + + // Sample reference points in WGS84 (not EPSG:3857) + const refPointsWgs84 = sampleReferencePointsInWgs84( + REF_POINTS_9, + tileCrsBounds, + this.projectTo4326, + ); + + // Project WGS84 points into globe common space via viewport.projectPosition + const refPointPositions: [number, number, number][] = []; + let csMinX = Number.POSITIVE_INFINITY; + let csMinY = Number.POSITIVE_INFINITY; + let csMaxX = Number.NEGATIVE_INFINITY; + let csMaxY = Number.NEGATIVE_INFINITY; + let latSum = 0; + + for (const [lng, lat] of refPointsWgs84) { + latSum += lat; + const posMin = project([lng, lat, minZ]); + refPointPositions.push([posMin[0]!, posMin[1]!, posMin[2]!]); + + if (posMin[0]! < csMinX) csMinX = posMin[0]!; + if (posMin[1]! < csMinY) csMinY = posMin[1]!; + if (posMin[0]! > csMaxX) csMaxX = posMin[0]!; + if (posMin[1]! > csMaxY) csMaxY = posMin[1]!; + + if (minZ !== maxZ) { + const posMax = project([lng, lat, maxZ]); + refPointPositions.push([posMax[0]!, posMax[1]!, posMax[2]!]); + + if (posMax[0]! < csMinX) csMinX = posMax[0]!; + if (posMax[1]! < csMinY) csMinY = posMax[1]!; + if (posMax[0]! > csMaxX) csMaxX = posMax[0]!; + if (posMax[1]! > csMaxY) csMaxY = posMax[1]!; + } + } + + const commonSpaceBounds: Bounds = [csMinX, csMinY, csMaxX, csMaxY]; + const centerLatitude = latSum / refPointsWgs84.length; + return { boundingVolume: makeOrientedBoundingBoxFromPoints(refPointPositions), commonSpaceBounds, + centerLatitude, }; } } @@ -534,6 +615,35 @@ function sampleReferencePointsInEPSG3857( return refPointPositions; } +/** + * Sample the selected reference points in WGS84 (EPSG:4326) + * + * Used for Globe view bounding volume computation where we need WGS84 + * coordinates instead of EPSG:3857. + * + * @param refPoints selected reference points. Each coordinate should be in [0-1] + * @param tileBounds the bounds of the tile in **tile CRS** [minX, minY, maxX, maxY] + * @param projectTo4326 projection function from tile CRS to WGS84 + */ +function sampleReferencePointsInWgs84( + refPoints: [number, number][], + tileBounds: [number, number, number, number], + projectTo4326: ProjectionFunction, +): [number, number][] { + const [minX, minY, maxX, maxY] = tileBounds; + const refPointPositions: [number, number][] = []; + + for (const [relX, relY] of refPoints) { + const geoX = minX + relX * (maxX - minX); + const geoY = minY + relY * (maxY - minY); + + const projected = projectTo4326(geoX, geoY); + refPointPositions.push(projected); + } + + return refPointPositions; +} + /** * Rescale positions from EPSG:3857 into deck.gl's common space * @@ -647,6 +757,7 @@ export function getTileIndices( maxZ: number; zRange: ZRange | null; projectTo3857: ProjectionFunction; + projectTo4326: ProjectionFunction; wgs84Bounds: CornerBounds; }, ): TileIndex[] { @@ -696,17 +807,33 @@ export function getTileIndices( // minZ to 0 const minZ = 0; - const { lowerLeft, upperRight } = wgs84Bounds; - const [minLng, minLat] = lowerLeft; - const [maxLng, maxLat] = upperRight; - const bottomLeft = lngLatToWorld([minLng, minLat]); - const topRight = lngLatToWorld([maxLng, maxLat]); - const bounds: Bounds = [ - bottomLeft[0], - bottomLeft[1], - topRight[0], - topRight[1], - ]; + // Convert WGS84 bounds to the appropriate common space for bounds filtering + let bounds: Bounds; + + if (project) { + // Globe view: project WGS84 bounds corners into globe common space + const { lowerLeft, upperRight } = wgs84Bounds; + const [minLng, minLat] = lowerLeft; + const [maxLng, maxLat] = upperRight; + const bl = project([minLng, minLat, 0]); + const tr = project([maxLng, maxLat, 0]); + const br = project([maxLng, minLat, 0]); + const tl = project([minLng, maxLat, 0]); + bounds = [ + Math.min(bl[0]!, tl[0]!, br[0]!, tr[0]!), + Math.min(bl[1]!, tl[1]!, br[1]!, tr[1]!), + Math.max(bl[0]!, tl[0]!, br[0]!, tr[0]!), + Math.max(bl[1]!, tl[1]!, br[1]!, tr[1]!), + ]; + } else { + // Mercator view: existing code + const { lowerLeft, upperRight } = wgs84Bounds; + const [minLng, minLat] = lowerLeft; + const [maxLng, maxLat] = upperRight; + const bottomLeft = lngLatToWorld([minLng, minLat]); + const topRight = lngLatToWorld([maxLng, maxLat]); + bounds = [bottomLeft[0], bottomLeft[1], topRight[0], topRight[1]]; + } // Start from coarsest overview const rootMatrix = metadata.tileMatrices[0]!; @@ -721,6 +848,7 @@ export function getTileIndices( new RasterTileNode(x, y, 0, { metadata, projectTo3857: opts.projectTo3857, + projectTo4326: opts.projectTo4326, }), ); } @@ -780,5 +908,11 @@ function getMetersPerPixelAtBoundingVolume( */ export const __TEST_EXPORTS = { computeProjectedTileBounds, + getOverlappingChildRange, + getMetersPerPixel, + getMetersPerPixelAtBoundingVolume, + rescaleEPSG3857ToCommonSpace, + sampleReferencePointsInEPSG3857, + sampleReferencePointsInWgs84, RasterTileNode, }; 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 8d3a3b61..66c7d85a 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 @@ -32,6 +32,7 @@ export class TileMatrixSetTileset extends Tileset2D { private tms: TileMatrixSet; private wgs84Bounds: CornerBounds; private projectTo3857: ProjectionFunction; + private projectTo4326: ProjectionFunction; constructor( opts: Tileset2DProps, @@ -47,6 +48,7 @@ export class TileMatrixSetTileset extends Tileset2D { super(opts); this.tms = tms; this.projectTo3857 = projectTo3857; + this.projectTo4326 = projectTo4326; if (!tms.boundingBox) { throw new Error( @@ -86,6 +88,7 @@ export class TileMatrixSetTileset extends Tileset2D { zRange: opts.zRange ?? null, wgs84Bounds: this.wgs84Bounds, projectTo3857: this.projectTo3857, + projectTo4326: this.projectTo4326, }); return tileIndices; diff --git a/packages/deck.gl-raster/tests/raster-tile-traversal.test.ts b/packages/deck.gl-raster/tests/raster-tile-traversal.test.ts new file mode 100644 index 00000000..5ace8aaa --- /dev/null +++ b/packages/deck.gl-raster/tests/raster-tile-traversal.test.ts @@ -0,0 +1,486 @@ +import { describe, expect, it } from "vitest"; +import _UTM31 from "../../morecantile/spec/schemas/tms/2.0/json/examples/tilematrixset/UTM31WGS84Quad.json"; +import _WebMercator from "../../morecantile/spec/schemas/tms/2.0/json/examples/tilematrixset/WebMercatorQuad.json"; +import type { + TileMatrix, + TileMatrixSet, +} from "../../morecantile/src/types/index"; +import { __TEST_EXPORTS } from "../src/raster-tileset/raster-tile-traversal"; +import type { ProjectionFunction } from "../src/raster-tileset/types"; + +const { + computeProjectedTileBounds, + getOverlappingChildRange, + getMetersPerPixel, + rescaleEPSG3857ToCommonSpace, + sampleReferencePointsInEPSG3857, + sampleReferencePointsInWgs84, + RasterTileNode, +} = __TEST_EXPORTS; + +const WebMercator = _WebMercator as TileMatrixSet; +const UTM31 = _UTM31 as TileMatrixSet; + +function findMatrix(tms: TileMatrixSet, id: string): TileMatrix { + const m = tms.tileMatrices.find((m) => m.id === id); + if (!m) throw new Error(`no matrix with id "${id}"`); + return m; +} + +// --------------------------------------------------------------------------- +// computeProjectedTileBounds +// --------------------------------------------------------------------------- +describe("computeProjectedTileBounds", () => { + it("returns correct bounds for WebMercatorQuad zoom 0 tile (0,0)", () => { + const matrix = findMatrix(WebMercator, "0"); + const bounds = computeProjectedTileBounds(matrix, { x: 0, y: 0 }); + // WebMercatorQuad zoom 0 has one tile covering the entire world + // EPSG:3857 full extent: ~[-20037508, -20037508, 20037508, 20037508] + const halfCirc = Math.PI * 6378137; + expect(bounds[0]).toBeCloseTo(-halfCirc, 0); + expect(bounds[1]).toBeCloseTo(-halfCirc, 0); + expect(bounds[2]).toBeCloseTo(halfCirc, 0); + expect(bounds[3]).toBeCloseTo(halfCirc, 0); + }); + + it("returns correct bounds for WebMercatorQuad zoom 1 tile (0,0)", () => { + const matrix = findMatrix(WebMercator, "1"); + const bounds = computeProjectedTileBounds(matrix, { x: 0, y: 0 }); + const halfCirc = Math.PI * 6378137; + // Top-left quadrant: [-halfCirc, 0, 0, halfCirc] + expect(bounds[0]).toBeCloseTo(-halfCirc, 0); + expect(bounds[1]).toBeCloseTo(0, 0); + expect(bounds[2]).toBeCloseTo(0, 0); + expect(bounds[3]).toBeCloseTo(halfCirc, 0); + }); + + it("returns correct bounds for UTM31 tile", () => { + // UTM31 matrix IDs start at "1", not "0" + const matrix = findMatrix(UTM31, "1"); + const bounds = computeProjectedTileBounds(matrix, { x: 0, y: 0 }); + // UTM31 should have bounds in meters, origin around (166021, ~9329005) + // Just verify it returns 4 finite numbers with min < max + expect(bounds).toHaveLength(4); + expect(bounds[0]).toBeLessThan(bounds[2]); // minX < maxX + expect(bounds[1]).toBeLessThan(bounds[3]); // minY < maxY + expect(Number.isFinite(bounds[0])).toBe(true); + }); +}); + +// --------------------------------------------------------------------------- +// rescaleEPSG3857ToCommonSpace +// --------------------------------------------------------------------------- +describe("rescaleEPSG3857ToCommonSpace", () => { + it("maps origin (0,0) in EPSG:3857 to center (256,256) in common space", () => { + const [x, y] = rescaleEPSG3857ToCommonSpace([0, 0]); + expect(x).toBeCloseTo(256, 5); + expect(y).toBeCloseTo(256, 5); + }); + + it("maps EPSG:3857 full extent to [0,512] range", () => { + const halfCirc = Math.PI * 6378137; + const [xMin, yMin] = rescaleEPSG3857ToCommonSpace([-halfCirc, -halfCirc]); + const [xMax, yMax] = rescaleEPSG3857ToCommonSpace([halfCirc, halfCirc]); + expect(xMin).toBeCloseTo(0, 5); + expect(yMin).toBeCloseTo(0, 5); + expect(xMax).toBeCloseTo(512, 5); + expect(yMax).toBeCloseTo(512, 5); + }); + + it("clamps Y values beyond Web Mercator bounds", () => { + const halfCirc = Math.PI * 6378137; + const beyondBounds = halfCirc * 2; + const [, yBeyond] = rescaleEPSG3857ToCommonSpace([0, beyondBounds]); + const [, yMax] = rescaleEPSG3857ToCommonSpace([0, halfCirc]); + // Should be clamped to the same value as halfCirc + expect(yBeyond).toBeCloseTo(yMax, 5); + }); +}); + +// --------------------------------------------------------------------------- +// sampleReferencePointsInEPSG3857 +// --------------------------------------------------------------------------- +describe("sampleReferencePointsInEPSG3857", () => { + it("identity projection returns input coordinates unchanged", () => { + const identity: ProjectionFunction = (x, y) => [x, y]; + const tileBounds: [number, number, number, number] = [100, 200, 300, 400]; + const refPoints: [number, number][] = [ + [0, 0], // lower-left corner + [1, 1], // upper-right corner + [0.5, 0.5], // center + ]; + const result = sampleReferencePointsInEPSG3857( + refPoints, + tileBounds, + identity, + ); + expect(result).toHaveLength(3); + // [0,0] → (100, 200) + expect(result[0]![0]).toBeCloseTo(100, 5); + expect(result[0]![1]).toBeCloseTo(200, 5); + // [1,1] → (300, 400) + expect(result[1]![0]).toBeCloseTo(300, 5); + expect(result[1]![1]).toBeCloseTo(400, 5); + // [0.5,0.5] → (200, 300) + expect(result[2]![0]).toBeCloseTo(200, 5); + expect(result[2]![1]).toBeCloseTo(300, 5); + }); +}); + +// --------------------------------------------------------------------------- +// getOverlappingChildRange +// --------------------------------------------------------------------------- +describe("getOverlappingChildRange", () => { + it("quadtree-like refinement: parent (0,0,z=0) covers 4 children", () => { + // WebMercatorQuad: z=0 is 1x1 tile, z=1 is 2x2 tiles + const parentMatrix = findMatrix(WebMercator, "0"); + const childMatrix = findMatrix(WebMercator, "1"); + const parentBounds = computeProjectedTileBounds(parentMatrix, { + x: 0, + y: 0, + }); + const range = getOverlappingChildRange(parentBounds, childMatrix); + expect(range.minCol).toBe(0); + expect(range.maxCol).toBe(1); + expect(range.minRow).toBe(0); + expect(range.maxRow).toBe(1); + }); + + it("quadtree-like refinement: z=1 tile (0,0) maps to z=2 quadrant", () => { + const parentMatrix = findMatrix(WebMercator, "1"); + const childMatrix = findMatrix(WebMercator, "2"); + const parentBounds = computeProjectedTileBounds(parentMatrix, { + x: 0, + y: 0, + }); + const range = getOverlappingChildRange(parentBounds, childMatrix); + expect(range.minCol).toBe(0); + // maxCol is 2 (not 1) because the parent boundary lands exactly on the + // child tile boundary, and Math.floor maps that to index 2 + expect(range.maxCol).toBe(2); + expect(range.minRow).toBe(0); + expect(range.maxRow).toBe(2); + }); + + it("z=1 tile (1,1) maps to z=2 lower-right quadrant", () => { + const parentMatrix = findMatrix(WebMercator, "1"); + const childMatrix = findMatrix(WebMercator, "2"); + const parentBounds = computeProjectedTileBounds(parentMatrix, { + x: 1, + y: 1, + }); + const range = getOverlappingChildRange(parentBounds, childMatrix); + expect(range.minCol).toBe(2); + expect(range.maxCol).toBe(3); + expect(range.minRow).toBe(2); + expect(range.maxRow).toBe(3); + }); +}); + +// --------------------------------------------------------------------------- +// getMetersPerPixel +// --------------------------------------------------------------------------- +describe("getMetersPerPixel", () => { + it("returns expected value at equator zoom 0", () => { + const earthCircumference = 40075016.686; + const expected = earthCircumference / 2 ** 8; // zoom 0, 2^(0+8) = 256 + const result = getMetersPerPixel(0, 0); + expect(result).toBeCloseTo(expected, 1); + }); + + it("decreases with increasing zoom", () => { + const z0 = getMetersPerPixel(0, 0); + const z1 = getMetersPerPixel(0, 1); + const z10 = getMetersPerPixel(0, 10); + expect(z0).toBeGreaterThan(z1); + expect(z1).toBeGreaterThan(z10); + // Each zoom level halves the meters per pixel + expect(z0 / z1).toBeCloseTo(2, 5); + }); + + it("decreases with increasing latitude (toward poles)", () => { + const equator = getMetersPerPixel(0, 5); + const lat60 = getMetersPerPixel(60, 5); + const lat80 = getMetersPerPixel(80, 5); + expect(equator).toBeGreaterThan(lat60); + expect(lat60).toBeGreaterThan(lat80); + // At 60° latitude, meters per pixel should be ~half of equator + expect(lat60 / equator).toBeCloseTo(0.5, 1); + }); +}); + +// --------------------------------------------------------------------------- +// RasterTileNode — insideBounds +// --------------------------------------------------------------------------- +describe("RasterTileNode.insideBounds", () => { + const identity: ProjectionFunction = (x, y) => [x, y]; + function makeNode( + x: number, + y: number, + z: number, + ): InstanceType { + return new RasterTileNode(x, y, z, { + metadata: WebMercator, + projectTo3857: identity, + projectTo4326: identity, + }); + } + + it("returns true for overlapping bounds", () => { + const node = makeNode(0, 0, 0); + const bounds = [0, 0, 300, 300] as [number, number, number, number]; + const commonSpaceBounds = [100, 100, 400, 400] as [ + number, + number, + number, + number, + ]; + expect(node.insideBounds(bounds, commonSpaceBounds)).toBe(true); + }); + + it("returns false for non-overlapping bounds", () => { + const node = makeNode(0, 0, 0); + const bounds = [0, 0, 50, 50] as [number, number, number, number]; + const commonSpaceBounds = [100, 100, 400, 400] as [ + number, + number, + number, + number, + ]; + expect(node.insideBounds(bounds, commonSpaceBounds)).toBe(false); + }); + + it("returns true for touching bounds (edge overlap)", () => { + const node = makeNode(0, 0, 0); + // Bounds touch at x=100: bounds goes up to 100, tile starts at 99 + const bounds = [0, 0, 100, 100] as [number, number, number, number]; + const commonSpaceBounds = [99, 0, 200, 200] as [ + number, + number, + number, + number, + ]; + expect(node.insideBounds(bounds, commonSpaceBounds)).toBe(true); + }); + + it("returns false for bounds that touch at exactly one edge (not overlapping)", () => { + const node = makeNode(0, 0, 0); + // tile starts exactly where bounds end — no overlap (< not <=) + const bounds = [0, 0, 100, 100] as [number, number, number, number]; + const commonSpaceBounds = [100, 0, 200, 200] as [ + number, + number, + number, + number, + ]; + expect(node.insideBounds(bounds, commonSpaceBounds)).toBe(false); + }); +}); + +// --------------------------------------------------------------------------- +// RasterTileNode — getBoundingVolume (Mercator path) +// --------------------------------------------------------------------------- +describe("RasterTileNode.getBoundingVolume (Mercator)", () => { + it("computes a bounding volume for WebMercatorQuad zoom 0", () => { + // Use identity projection (pretend tile CRS is already EPSG:3857) + const identity: ProjectionFunction = (x, y) => [x, y]; + const node = new RasterTileNode(0, 0, 0, { + metadata: WebMercator, + projectTo3857: identity, + projectTo4326: identity, + }); + + const zRange: [number, number] = [0, 0]; + const { boundingVolume, commonSpaceBounds } = node.getBoundingVolume( + zRange, + null, + ); + + // Should have a valid OrientedBoundingBox + expect(boundingVolume).toBeDefined(); + expect(boundingVolume.center).toBeDefined(); + expect(boundingVolume.halfAxes).toBeDefined(); + + // Common space bounds should span most of [0, 512] + const [minX, minY, maxX, maxY] = commonSpaceBounds; + expect(maxX - minX).toBeGreaterThan(400); // Should be ~512 wide + expect(maxY - minY).toBeGreaterThan(400); // Should be ~512 tall + }); + + it("z=1 tiles have smaller bounding volumes than z=0", () => { + const identity: ProjectionFunction = (x, y) => [x, y]; + + const nodeZ0 = new RasterTileNode(0, 0, 0, { + metadata: WebMercator, + projectTo3857: identity, + projectTo4326: identity, + }); + const nodeZ1 = new RasterTileNode(0, 0, 1, { + metadata: WebMercator, + projectTo3857: identity, + projectTo4326: identity, + }); + + const zRange: [number, number] = [0, 0]; + const { commonSpaceBounds: csZ0 } = nodeZ0.getBoundingVolume(zRange, null); + const { commonSpaceBounds: csZ1 } = nodeZ1.getBoundingVolume(zRange, null); + + const widthZ0 = csZ0[2] - csZ0[0]; + const widthZ1 = csZ1[2] - csZ1[0]; + expect(widthZ0).toBeGreaterThan(widthZ1); + }); +}); + +// --------------------------------------------------------------------------- +// RasterTileNode — children +// --------------------------------------------------------------------------- +describe("RasterTileNode.children", () => { + it("WebMercatorQuad z=0 tile has 4 children at z=1", () => { + const identity: ProjectionFunction = (x, y) => [x, y]; + const node = new RasterTileNode(0, 0, 0, { + metadata: WebMercator, + projectTo3857: identity, + projectTo4326: identity, + }); + + const children = node.children; + expect(children).not.toBeNull(); + expect(children).toHaveLength(4); + + // Children should be at z=1 + for (const child of children!) { + expect(child.z).toBe(1); + } + + // Should cover all 4 quadrants + const coords = children!.map((c) => [c.x, c.y]); + expect(coords).toContainEqual([0, 0]); + expect(coords).toContainEqual([1, 0]); + expect(coords).toContainEqual([0, 1]); + expect(coords).toContainEqual([1, 1]); + }); + + it("finest zoom level has no children", () => { + const identity: ProjectionFunction = (x, y) => [x, y]; + const maxZ = WebMercator.tileMatrices.length - 1; + const node = new RasterTileNode(0, 0, maxZ, { + metadata: WebMercator, + projectTo3857: identity, + projectTo4326: identity, + }); + + expect(node.children).toBeNull(); + }); +}); + +// --------------------------------------------------------------------------- +// sampleReferencePointsInWgs84 +// --------------------------------------------------------------------------- +describe("sampleReferencePointsInWgs84", () => { + it("identity projection returns input coordinates unchanged", () => { + const identity: ProjectionFunction = (x, y) => [x, y]; + const tileBounds: [number, number, number, number] = [100, 200, 300, 400]; + const refPoints: [number, number][] = [ + [0, 0], + [1, 1], + [0.5, 0.5], + ]; + const result = sampleReferencePointsInWgs84( + refPoints, + tileBounds, + identity, + ); + expect(result).toHaveLength(3); + expect(result[0]![0]).toBeCloseTo(100, 5); + expect(result[0]![1]).toBeCloseTo(200, 5); + expect(result[1]![0]).toBeCloseTo(300, 5); + expect(result[1]![1]).toBeCloseTo(400, 5); + expect(result[2]![0]).toBeCloseTo(200, 5); + expect(result[2]![1]).toBeCloseTo(300, 5); + }); +}); + +// --------------------------------------------------------------------------- +// RasterTileNode — getBoundingVolume (Globe path) +// --------------------------------------------------------------------------- +describe("RasterTileNode.getBoundingVolume (Globe)", () => { + it("computes a bounding volume using the project function", () => { + const identity: ProjectionFunction = (x, y) => [x, y]; + + // Mock globe project function that maps [lng, lat, z] to 3D common space + // Simple sphere: x = cos(lat)*cos(lng), y = cos(lat)*sin(lng), z = sin(lat) + // But for testing, a simple linear transform is sufficient to verify + // the plumbing works + const mockProject = (xyz: number[]): number[] => { + return [xyz[0]! * 10, xyz[1]! * 10, xyz[2]! || 0]; + }; + + const node = new RasterTileNode(0, 0, 0, { + metadata: WebMercator, + projectTo3857: identity, + projectTo4326: identity, + }); + + const zRange: [number, number] = [0, 0]; + const { boundingVolume, commonSpaceBounds, centerLatitude } = + node.getBoundingVolume(zRange, mockProject); + + // Should have a valid OrientedBoundingBox + expect(boundingVolume).toBeDefined(); + expect(boundingVolume.center).toBeDefined(); + expect(boundingVolume.halfAxes).toBeDefined(); + + // Common space bounds should be defined + const [minX, minY, maxX, maxY] = commonSpaceBounds; + expect(Number.isFinite(minX)).toBe(true); + expect(Number.isFinite(minY)).toBe(true); + expect(maxX).toBeGreaterThan(minX); + expect(maxY).toBeGreaterThan(minY); + + // Center latitude should be finite + expect(Number.isFinite(centerLatitude)).toBe(true); + }); + + it("produces different bounding volumes than the Mercator path", () => { + const identity: ProjectionFunction = (x, y) => [x, y]; + const mockProject = (xyz: number[]): number[] => { + return [xyz[0]! * 5, xyz[1]! * 5, 0]; + }; + + const node = new RasterTileNode(0, 0, 0, { + metadata: WebMercator, + projectTo3857: identity, + projectTo4326: identity, + }); + + const zRange: [number, number] = [0, 0]; + const mercator = node.getBoundingVolume(zRange, null); + const globe = node.getBoundingVolume(zRange, mockProject); + + // The bounding volumes should differ because the common spaces differ + expect(globe.commonSpaceBounds[0]).not.toBeCloseTo( + mercator.commonSpaceBounds[0], + 1, + ); + }); + + it("returns centerLatitude from WGS84 reference points", () => { + // projectTo4326 that always returns lng=0, lat=45 + const mockTo4326: ProjectionFunction = (_x, _y) => [0, 45]; + const identity: ProjectionFunction = (x, y) => [x, y]; + const mockProject = (xyz: number[]): number[] => [ + xyz[0]!, + xyz[1]!, + xyz[2]! || 0, + ]; + + const node = new RasterTileNode(0, 0, 0, { + metadata: WebMercator, + projectTo3857: identity, + projectTo4326: mockTo4326, + }); + + const { centerLatitude } = node.getBoundingVolume([0, 0], mockProject); + expect(centerLatitude).toBeCloseTo(45, 5); + }); +}); 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); + } + }); +}); diff --git a/pnpm-lock.yaml b/pnpm-lock.yaml index c62c2090..231685c6 100644 --- a/pnpm-lock.yaml +++ b/pnpm-lock.yaml @@ -102,6 +102,64 @@ importers: specifier: ^7.3.1 version: 7.3.1(@types/node@25.1.0)(tsx@4.21.0) + examples/globe-view: + dependencies: + '@deck.gl/core': + specifier: ^9.2.8 + version: 9.2.8 + '@deck.gl/geo-layers': + specifier: ^9.2.8 + version: 9.2.8(@deck.gl/core@9.2.8)(@deck.gl/extensions@9.2.5(@deck.gl/core@9.2.8)(@luma.gl/core@9.2.6)(@luma.gl/engine@9.2.6(@luma.gl/core@9.2.6)(@luma.gl/shadertools@9.2.6(@luma.gl/core@9.2.6))))(@deck.gl/layers@9.2.8(@deck.gl/core@9.2.8)(@loaders.gl/core@4.3.4)(@luma.gl/core@9.2.6)(@luma.gl/engine@9.2.6(@luma.gl/core@9.2.6)(@luma.gl/shadertools@9.2.6(@luma.gl/core@9.2.6))))(@deck.gl/mesh-layers@9.2.8(@deck.gl/core@9.2.8)(@loaders.gl/core@4.3.4)(@luma.gl/core@9.2.6)(@luma.gl/engine@9.2.6(@luma.gl/core@9.2.6)(@luma.gl/shadertools@9.2.6(@luma.gl/core@9.2.6)))(@luma.gl/gltf@9.2.6(@luma.gl/constants@9.2.6)(@luma.gl/core@9.2.6)(@luma.gl/engine@9.2.6(@luma.gl/core@9.2.6)(@luma.gl/shadertools@9.2.6(@luma.gl/core@9.2.6)))(@luma.gl/shadertools@9.2.6(@luma.gl/core@9.2.6)))(@luma.gl/shadertools@9.2.6(@luma.gl/core@9.2.6)))(@loaders.gl/core@4.3.4)(@luma.gl/constants@9.2.6)(@luma.gl/core@9.2.6)(@luma.gl/engine@9.2.6(@luma.gl/core@9.2.6)(@luma.gl/shadertools@9.2.6(@luma.gl/core@9.2.6))) + '@deck.gl/layers': + specifier: ^9.2.8 + version: 9.2.8(@deck.gl/core@9.2.8)(@loaders.gl/core@4.3.4)(@luma.gl/core@9.2.6)(@luma.gl/engine@9.2.6(@luma.gl/core@9.2.6)(@luma.gl/shadertools@9.2.6(@luma.gl/core@9.2.6))) + '@deck.gl/mesh-layers': + specifier: ^9.2.8 + version: 9.2.8(@deck.gl/core@9.2.8)(@loaders.gl/core@4.3.4)(@luma.gl/core@9.2.6)(@luma.gl/engine@9.2.6(@luma.gl/core@9.2.6)(@luma.gl/shadertools@9.2.6(@luma.gl/core@9.2.6)))(@luma.gl/gltf@9.2.6(@luma.gl/constants@9.2.6)(@luma.gl/core@9.2.6)(@luma.gl/engine@9.2.6(@luma.gl/core@9.2.6)(@luma.gl/shadertools@9.2.6(@luma.gl/core@9.2.6)))(@luma.gl/shadertools@9.2.6(@luma.gl/core@9.2.6)))(@luma.gl/shadertools@9.2.6(@luma.gl/core@9.2.6)) + '@deck.gl/react': + specifier: ^9.2.7 + version: 9.2.9(@deck.gl/core@9.2.8)(@deck.gl/widgets@9.2.9(@deck.gl/core@9.2.8)(@luma.gl/core@9.2.6))(react-dom@19.2.4(react@19.2.4))(react@19.2.4) + '@developmentseed/deck.gl-geotiff': + specifier: workspace:^ + version: link:../../packages/deck.gl-geotiff + '@developmentseed/deck.gl-raster': + specifier: workspace:^ + version: link:../../packages/deck.gl-raster + '@developmentseed/geotiff': + specifier: workspace:^ + version: link:../../packages/geotiff + '@luma.gl/core': + specifier: ^9.2.6 + version: 9.2.6 + '@luma.gl/shadertools': + specifier: ^9.2.6 + version: 9.2.6(@luma.gl/core@9.2.6) + '@luma.gl/webgl': + specifier: ^9.2.6 + version: 9.2.6(@luma.gl/core@9.2.6) + proj4: + specifier: ^2.20.2 + version: 2.20.2 + react: + specifier: ^19.2.4 + version: 19.2.4 + react-dom: + specifier: ^19.2.4 + version: 19.2.4(react@19.2.4) + devDependencies: + '@types/react': + specifier: ^19.2.10 + version: 19.2.10 + '@types/react-dom': + specifier: ^19.2.3 + version: 19.2.3(@types/react@19.2.10) + '@vitejs/plugin-react': + specifier: ^5.1.2 + version: 5.1.2(vite@7.3.1(@types/node@25.1.0)(tsx@4.21.0)) + vite: + specifier: ^7.3.1 + version: 7.3.1(@types/node@25.1.0)(tsx@4.21.0) + examples/land-cover: dependencies: '@deck.gl/core': @@ -734,6 +792,20 @@ packages: '@luma.gl/gltf': ^9.2.6 '@luma.gl/shadertools': ^9.2.6 + '@deck.gl/react@9.2.9': + resolution: {integrity: sha512-ZADiRJhT8dI1z6NfC3cJ62j8nt4j/5AfUmXSHNn6hnOTJFJB1YScERnmM2lz9LYBwjlWTWEeTDIxQ66S4nClYg==} + peerDependencies: + '@deck.gl/core': ^9.2.8 + '@deck.gl/widgets': ^9.2.8 + react: '>=16.3.0' + react-dom: '>=16.3.0' + + '@deck.gl/widgets@9.2.9': + resolution: {integrity: sha512-dFIT1sJZ8gxZE7l+b62TY5LH/92ABcsukK7jH3L3T6PYlbFOnNCby6dkZRySYbhWPE77lykXXsx1V/NQQlqs8Q==} + peerDependencies: + '@deck.gl/core': ^9.2.8 + '@luma.gl/core': ^9.2.6 + '@developmentseed/lzw-tiff-decoder@0.2.2': resolution: {integrity: sha512-bsBIdV1LyqcrtnYtIYu3C/297X2wxeYkmmhHzR02OOX823sGxT8GLGdenTmwoXvNWkozk3b+ptXtJeSUmNPaTA==} @@ -2178,6 +2250,9 @@ packages: potpack@2.1.0: resolution: {integrity: sha512-pcaShQc1Shq0y+E7GqJqvZj8DTthWV1KeHGdi0Z6IAin2Oi3JnLCOfwnCo84qc+HAp52wT9nK9H7FAJp5a44GQ==} + preact@10.28.4: + resolution: {integrity: sha512-uKFfOHWuSNpRFVTnljsCluEFq57OKT+0QdOiQo8XWnQ/pSvg7OpX5eNOejELXJMWy+BwM2nobz0FkvzmnpCNsQ==} + prettier@3.8.1: resolution: {integrity: sha512-UOnG6LftzbdaHZcKoPFtOcCKztrQ57WkHDeRD9t/PTQtmT0NHSeWWepj6pS0z/N7+08BHFDQVUrfmfMRcZwbMg==} engines: {node: '>=14'} @@ -2955,6 +3030,19 @@ snapshots: transitivePeerDependencies: - '@loaders.gl/core' + '@deck.gl/react@9.2.9(@deck.gl/core@9.2.8)(@deck.gl/widgets@9.2.9(@deck.gl/core@9.2.8)(@luma.gl/core@9.2.6))(react-dom@19.2.4(react@19.2.4))(react@19.2.4)': + dependencies: + '@deck.gl/core': 9.2.8 + '@deck.gl/widgets': 9.2.9(@deck.gl/core@9.2.8)(@luma.gl/core@9.2.6) + react: 19.2.4 + react-dom: 19.2.4(react@19.2.4) + + '@deck.gl/widgets@9.2.9(@deck.gl/core@9.2.8)(@luma.gl/core@9.2.6)': + dependencies: + '@deck.gl/core': 9.2.8 + '@luma.gl/core': 9.2.6 + preact: 10.28.4 + '@developmentseed/lzw-tiff-decoder@0.2.2': {} '@esbuild/aix-ppc64@0.27.2': @@ -4342,6 +4430,8 @@ snapshots: potpack@2.1.0: {} + preact@10.28.4: {} + prettier@3.8.1: {} process-nextick-args@2.0.1: {}