From 8d8c98ebb1b92fc40233c4c039bdbde4e261a96a Mon Sep 17 00:00:00 2001 From: Kyle Barron Date: Mon, 15 Dec 2025 19:56:25 -0500 Subject: [PATCH 01/14] Add in existing code --- .../raster-tileset/raster-tile-traversal.ts | 535 ++++++++++++++++++ .../src/raster-tileset/raster-tileset-2d.ts | 344 +++++++++++ 2 files changed, 879 insertions(+) create mode 100644 packages/deck.gl-raster/src/raster-tileset/raster-tile-traversal.ts create mode 100644 packages/deck.gl-raster/src/raster-tileset/raster-tileset-2d.ts 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 new file mode 100644 index 0000000..c6f17fb --- /dev/null +++ b/packages/deck.gl-raster/src/raster-tileset/raster-tile-traversal.ts @@ -0,0 +1,535 @@ +/** + * This file implements tile traversal for generic 2D tilesets defined by COG + * tile layouts. + * + * The main algorithm works as follows: + * 1. Start at the root tile(s) (z=0, covers the entire image, but not + * necessarily the whole world) + * 2. Test if each tile is visible using viewport frustum culling + * 3. For visible tiles, compute distance-based LOD (Level of Detail) + * 4. If LOD is insufficient, recursively subdivide into 4 child tiles + * 5. Select tiles at appropriate zoom levels based on distance from camera + * + * The result is a set of tiles at varying zoom levels that efficiently + * cover the visible area with appropriate detail. + */ + +import { + _GlobeViewport, + assert, + Viewport, + WebMercatorViewport, +} from "@deck.gl/core"; +import { + CullingVolume, + Plane, + makeOrientedBoundingBoxFromPoints, +} from "@math.gl/culling"; + +import type { + TileMatrixSet, + TileMatrix, + TileIndex, + ZRange, +} from "../raster-tileset/types.js"; + +/** + * The size of the entire world in deck.gl's common coordinate space. + * + * The world always spans [0, 512] in both X and Y in Web Mercator common space. + * + * The origin (0,0) is at the top-left corner, and (512,512) is at the + * bottom-right. + */ +const WORLD_SIZE = 512; + +// Reference points used to sample tile boundaries for bounding volume +// calculation. +// +// In upstream deck.gl code, such reference points are only used in non-Web +// Mercator projections because the OSM tiling scheme is designed for Web +// Mercator and the OSM tile extents are already in Web Mercator projection. So +// using Axis-Aligned bounding boxes based on tile extents is sufficient for +// frustum culling in Web Mercator viewports. +// +// In upstream code these reference points are used for Globe View where the OSM +// tile indices _projected into longitude-latitude bounds in Globe View space_ +// are no longer axis-aligned, and oriented bounding boxes must be used instead. +// +// In the context of generic tiling grids which are often not in Web Mercator +// projection, we must use the reference points approach because the grid tiles +// will never be exact axis aligned boxes in Web Mercator space. + +// For most tiles: sample 4 corners and center (5 points total) +const REF_POINTS_5: [number, number][] = [ + [0.5, 0.5], // center + [0, 0], // top-left + [0, 1], // bottom-left + [1, 0], // top-right + [1, 1], // bottom-right +]; + +// For higher detail: add 4 edge midpoints (9 points total) +const REF_POINTS_9 = REF_POINTS_5.concat([ + [0, 0.5], // left edge + [0.5, 0], // top edge + [1, 0.5], // right edge + [0.5, 1], // bottom edge +]); + +/** + * Raster Tile Node - similar to OSMNode but for a generic raster tileset's + * tile structure. + * + * Represents a single tile in the COG internal tiling pyramid. + * + * COG tile nodes use the following coordinate system: + * + * - x: tile column (0 to TileMatrix.tilesX, left to right) + * - y: tile row (0 to TileMatrix.tilesY, top to bottom) + * - z: overview level. This uses TileMatrixSet ordering where: 0 = coarsest, higher = finer + */ +export class RasterTileNode { + /** Index across a row */ + x: number; + /** Index down a column */ + y: number; + /** TileMatrixSet-style zoom index (higher = finer detail) */ + z: number; + + private metadata: TileMatrixSet; + + /** + * Flag indicating whether any descendant of this tile is visible. + * + * Used to prevent loading parent tiles when children are visible (avoids + * overdraw). + */ + private childVisible?: boolean; + + /** + * Flag indicating this tile should be rendered + * + * Set to `true` when this is the appropriate LOD for its distance from camera. + */ + private selected?: boolean; + + /** A cache of the children of this node. */ + private _children?: RasterTileNode[]; + + constructor(x: number, y: number, z: number, metadata: TileMatrixSet) { + this.x = x; + this.y = y; + this.z = z; + this.metadata = metadata; + } + + /** Get overview info for this tile's z level */ + get overview(): TileMatrix { + return this.metadata.tileMatrices[this.z]!; + } + + /** Get the children of this node. */ + get children(): RasterTileNode[] { + if (!this._children) { + const maxZ = this.metadata.tileMatrices.length - 1; + if (this.z >= maxZ) { + // Already at finest resolution, no children + return []; + } + + // In TileMatrixSet ordering: refine to z + 1 (finer detail) + const childZ = this.z + 1; + const parentOverview = this.overview; + const childOverview = this.metadata.tileMatrices[childZ]!; + + // Calculate scale factor between levels + const scaleFactor = parentOverview.cellSize / childOverview.cellSize; + + // Generate child tiles + this._children = []; + for (let dy = 0; dy < scaleFactor; dy++) { + for (let dx = 0; dx < scaleFactor; dx++) { + const childX = this.x * scaleFactor + dx; + const childY = this.y * scaleFactor + dy; + + // Only create child if it's within bounds + // Some tiles on the edges might not need to be created at higher + // resolutions (higher map zoom level) + if ( + childX < childOverview.tileWidth && + childY < childOverview.tileHeight + ) { + this._children.push( + new RasterTileNode(childX, childY, childZ, this.metadata), + ); + } + } + } + } + return this._children; + } + + /** + * Update tile visibility using frustum culling + * This follows the pattern from OSMNode + */ + update(params: { + viewport: Viewport; + project: ((xyz: number[]) => number[]) | null; + cullingVolume: CullingVolume; + elevationBounds: ZRange; + /** Minimum (coarsest) COG overview level */ + minZ: number; + /** Maximum (finest) COG overview level */ + maxZ?: number; + }): boolean { + const { + viewport, + cullingVolume, + elevationBounds, + minZ, + maxZ = this.metadata.tileMatrices.length - 1, + project, + } = params; + + // Get bounding volume for this tile + const boundingVolume = this.getBoundingVolume(elevationBounds, project); + + // Note: this is a part of the upstream code because they have _generic_ + // tiling systems, where the client doesn't know whether a given xyz tile + // actually exists. So the idea of `bounds` is to avoid even trying to fetch + // tiles that the user doesn't care about (think oceans) + // + // But in our case, we have known bounds from the COG metadata. So the tiles + // are explicitly constructed to match only tiles that exist. + + // Check if tile is within user-specified bounds + // if (bounds && !this.insideBounds(bounds)) { + // return false; + // } + + console.log("=== FRUSTUM CULLING DEBUG ==="); + console.log(`Tile: ${this.x}, ${this.y}, ${this.z}`); + console.log("Bounding volume center:", boundingVolume.center); + console.log("Bounding volume halfSize:", boundingVolume.halfSize); + console.log("Viewport cameraPosition:", viewport.cameraPosition); + console.log( + "Viewport pitch:", + viewport instanceof WebMercatorViewport ? viewport.pitch : "N/A", + ); + + for (let i = 0; i < cullingVolume.planes.length; i++) { + const plane = cullingVolume.planes[i]!; + const result = boundingVolume.intersectPlane(plane); + const planeNames = ["left", "right", "bottom", "top", "near", "far"]; + + // Calculate signed distance from OBB center to plane + const centerDist = + plane.normal.x * boundingVolume.center.x + + plane.normal.y * boundingVolume.center.y + + plane.normal.z * boundingVolume.center.z + + plane.distance; + + console.log( + `Plane ${i} (${planeNames[i]}): normal=[${plane.normal.x.toFixed(3)}, ${plane.normal.y.toFixed(3)}, ${plane.normal.z.toFixed(3)}], ` + + `distance=${plane.distance.toFixed(3)}, centerDist=${centerDist.toFixed(3)}, result=${result} (${result === 1 ? "INSIDE" : result === 0 ? "INTERSECT" : "OUTSIDE"})`, + ); + } + console.log("=== END FRUSTUM DEBUG ==="); + + // Frustum culling + // Test if tile's bounding volume intersects the camera frustum + // Returns: <0 if outside, 0 if intersecting, >0 if fully inside + const isInside = cullingVolume.computeVisibility(boundingVolume); + console.log( + `Tile ${this.x},${this.y},${this.z} frustum check: ${isInside} (${isInside < 0 ? "CULLED" : "VISIBLE"})`, + ); + if (isInside < 0) { + return false; + } + + // LOD (Level of Detail) selection + // Only select this tile if no child is visible (prevents overlapping tiles) + if (!this.childVisible) { + let { z } = this; + + if (z < maxZ && z >= minZ) { + // Compute distance-based LOD adjustment + // Tiles farther from camera can use lower zoom levels (larger tiles) + // Distance is normalized by viewport height to be resolution-independent + const distance = + (boundingVolume.distanceTo(viewport.cameraPosition) * + viewport.scale) / + viewport.height; + // Increase effective zoom level based on log2(distance) + // e.g., if distance=4, accept tiles 2 levels lower than maxZ + z += Math.floor(Math.log2(distance)); + } + + if (z >= maxZ) { + // This tile's LOD is sufficient for its distance - select it for rendering + this.selected = true; + return true; + } + } + + // LOD is not enough, recursively test child tiles + this.selected = false; + this.childVisible = true; + + for (const child of this.children) { + child.update(params); + } + + // // NOTE: this deviates from upstream; we could move to the upstream code if + // // we pass in maxZ correctly I think + // if (children.length === 0) { + // // No children available (at finest resolution), select this tile + // this.selected = true; + // return true; + // } + + // for (const child of children) { + // child.update(params); + // } + return true; + } + + /** + * Collect all tiles marked as selected in the tree. + * Recursively traverses the entire tree and gathers tiles where selected=true. + * + * @param result - Accumulator array for selected tiles + * @returns Array of selected OSMNode tiles + */ + getSelected(result: RasterTileNode[] = []): RasterTileNode[] { + if (this.selected) { + result.push(this); + } + if (this._children) { + for (const node of this._children) { + node.getSelected(result); + } + } + return result; + } + + /** + * Calculate the 3D bounding volume for this tile in deck.gl's common + * coordinate space for frustum culling. + * + */ + getBoundingVolume( + zRange: ZRange, + project: ((xyz: number[]) => number[]) | null, + ) { + const overview = this.overview; + const { tileMatrices } = this.metadata; + const { tileWidth, tileHeight } = tileMatrices[tileMatrices.length - 1]!; + + // Use geotransform to calculate tile bounds + // geotransform: [a, b, c, d, e, f] where: + // x_geo = a * col + b * row + c + // y_geo = d * col + e * row + f + const [a, b, c, d, e, f] = overview.geotransform; + + // Calculate pixel coordinates for this tile's extent + const pixelMinCol = this.x * tileWidth; + const pixelMinRow = this.y * tileHeight; + const pixelMaxCol = (this.x + 1) * tileWidth; + const pixelMaxRow = (this.y + 1) * tileHeight; + + // Sample reference points across the tile surface + const refPoints = REF_POINTS_9; + + console.log("refPoints", refPoints); + + /** Reference points positions in image CRS */ + const refPointPositionsImage: [number, number][] = []; + + for (const [pX, pY] of refPoints) { + // pX, pY are in [0, 1] range + // Interpolate pixel coordinates within the tile + const col = pixelMinCol + pX * (pixelMaxCol - pixelMinCol); + const row = pixelMinRow + pY * (pixelMaxRow - pixelMinRow); + + // Convert pixel coordinates to geographic coordinates using geotransform + const geoX = a * col + b * row + c; + const geoY = d * col + e * row + f; + + refPointPositionsImage.push([geoX, geoY]); + } + + console.log("refPointPositionsImage (image CRS):", refPointPositionsImage); + console.log("Geotransform [a,b,c,d,e,f]:", [a, b, c, d, e, f]); + + if (project) { + assert( + false, + "TODO: implement bounding volume implementation in Globe view", + ); + // Reproject positions to wgs84 instead, then pass them into `project` + // return makeOrientedBoundingBoxFromPoints(refPointPositions); + } + + /** Reference points positions in EPSG 3857 */ + const refPointPositionsProjected: [number, number][] = []; + + for (const [pX, pY] of refPointPositionsImage) { + // Reproject to Web Mercator (EPSG 3857) + const projected = this.metadata.projectTo3857([pX, pY]); + refPointPositionsProjected.push(projected); + + // Also log WGS84 for comparison + const wgs84 = this.metadata.projectToWgs84([pX, pY]); + console.log( + `Image [${pX.toFixed(2)}, ${pY.toFixed(2)}] -> WGS84 [${wgs84[0].toFixed(6)}, ${wgs84[1].toFixed(6)}] -> WebMerc [${projected[0].toFixed(2)}, ${projected[1].toFixed(2)}]`, + ); + } + + console.log( + "refPointPositionsProjected (EPSG:3857):", + refPointPositionsProjected, + ); + + // Convert from Web Mercator meters to deck.gl's common space (world units) + // Web Mercator range: [-20037508.34, 20037508.34] meters + // deck.gl world space: [0, 512] + const WEB_MERCATOR_MAX = 20037508.342789244; // Half Earth circumference + + /** Reference points positions in deck.gl world space */ + const refPointPositionsWorld: [number, number, number][] = []; + + for (const [mercX, mercY] of refPointPositionsProjected) { + // X: offset from [-20M, 20M] to [0, 40M], then normalize to [0, 512] + const worldX = + ((mercX + WEB_MERCATOR_MAX) / (2 * WEB_MERCATOR_MAX)) * WORLD_SIZE; + + // Y: same transformation WITHOUT flip + // Testing hypothesis: Y-flip might be incorrect since geotransform already handles orientation + const worldY = + ((mercY + WEB_MERCATOR_MAX) / (2 * WEB_MERCATOR_MAX)) * WORLD_SIZE; + + console.log( + `WebMerc [${mercX.toFixed(2)}, ${mercY.toFixed(2)}] -> World [${worldX.toFixed(4)}, ${worldY.toFixed(4)}]`, + ); + + // Add z-range minimum + refPointPositionsWorld.push([worldX, worldY, zRange[0]]); + } + + // Add top z-range if elevation varies + if (zRange[0] !== zRange[1]) { + for (const [mercX, mercY] of refPointPositionsProjected) { + const worldX = + ((mercX + WEB_MERCATOR_MAX) / (2 * WEB_MERCATOR_MAX)) * WORLD_SIZE; + const worldY = + WORLD_SIZE - + ((mercY + WEB_MERCATOR_MAX) / (2 * WEB_MERCATOR_MAX)) * WORLD_SIZE; + + refPointPositionsWorld.push([worldX, worldY, zRange[1]]); + } + } + + console.log("refPointPositionsWorld", refPointPositionsWorld); + console.log("zRange used:", zRange); + + const obb = makeOrientedBoundingBoxFromPoints(refPointPositionsWorld); + console.log("Created OBB center:", obb.center); + console.log("Created OBB halfAxes:", obb.halfAxes); + + return obb; + } + + /** + * Convert COG coordinates to lng/lat + * This is a placeholder - needs proper projection library (proj4js) + */ + private cogCoordsToLngLat([x, y]: [number, number]): number[] { + const [lng, lat] = this.metadata.projectToWgs84([x, y]); + return [lng, lat, 0]; + } +} + +/** + * Get tile indices visible in viewport + * Uses frustum culling similar to OSM implementation + * + * Overviews follow TileMatrixSet ordering: index 0 = coarsest, higher = finer + */ +export function getTileIndices( + metadata: TileMatrixSet, + opts: { + viewport: Viewport; + maxZ: number; + zRange: ZRange | null; + }, +): TileIndex[] { + const { viewport, maxZ, zRange } = opts; + + // console.log("=== getTileIndices called ==="); + // console.log("Viewport:", viewport); + // console.log("maxZ:", maxZ); + // console.log("COG metadata tileMatrices count:", metadata.tileMatrices.length); + // console.log("COG bbox:", metadata.bbox); + + const project: ((xyz: number[]) => number[]) | null = + viewport instanceof _GlobeViewport && viewport.resolution + ? viewport.projectPosition + : null; + + // Get the culling volume of the current camera + const planes: Plane[] = Object.values(viewport.getFrustumPlanes()).map( + ({ normal, distance }) => new Plane(normal.clone().negate(), distance), + ); + const cullingVolume = new CullingVolume(planes); + + // Project zRange from meters to common space + const unitsPerMeter = viewport.distanceScales.unitsPerMeter[2]!; + const elevationMin = (zRange && zRange[0] * unitsPerMeter) || 0; + const elevationMax = (zRange && zRange[1] * unitsPerMeter) || 0; + + // Optimization: For low-pitch views, only consider tiles at maxZ level + // At low pitch (top-down view), all tiles are roughly the same distance, + // so we don't need the LOD pyramid - just use the finest level + const minZ = + viewport instanceof WebMercatorViewport && viewport.pitch <= 60 ? maxZ : 0; + + // Start from coarsest overview + const coarsestOverview = metadata.tileMatrices[0]!; + + // Create root tiles at coarsest level + // In contrary to OSM tiling, we might have more than one tile at the + // coarsest level (z=0) + const roots: RasterTileNode[] = []; + for (let y = 0; y < coarsestOverview.tileHeight; y++) { + for (let x = 0; x < coarsestOverview.tileWidth; x++) { + roots.push(new RasterTileNode(x, y, 0, metadata)); + } + } + + // Traverse and update visibility + const traversalParams = { + viewport, + project, + cullingVolume, + elevationBounds: [elevationMin, elevationMax] as ZRange, + minZ, + maxZ, + }; + console.log("Traversal params:", traversalParams); + + for (const root of roots) { + root.update(traversalParams); + } + console.log("roots", roots); + + // Collect selected tiles + const selectedNodes: RasterTileNode[] = []; + for (const root of roots) { + root.getSelected(selectedNodes); + } + + return selectedNodes; +} 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 new file mode 100644 index 0000000..e3e8ccf --- /dev/null +++ b/packages/deck.gl-raster/src/raster-tileset/raster-tileset-2d.ts @@ -0,0 +1,344 @@ +/** + * COGTileset2D - Improved Implementation with Frustum Culling + * + * This version properly implements frustum culling and bounding volume calculations + * following the pattern from deck.gl's OSM tile indexing. + */ + +import { Viewport, WebMercatorViewport } from "@deck.gl/core"; +import { _Tileset2D as Tileset2D } from "@deck.gl/geo-layers"; +import type { Tileset2DProps } from "@deck.gl/geo-layers/dist/tileset-2d"; +import type { ZRange } from "@deck.gl/geo-layers/dist/tileset-2d/types"; +import { Matrix4 } from "@math.gl/core"; +import { GeoTIFF, GeoTIFFImage } from "geotiff"; +import proj4 from "proj4"; + +import { getTileIndices } from "./raster-tile-traversal"; +import type { COGMetadata, COGTileIndex, COGOverview, Bounds } from "./types"; + +const OGC_84 = { + $schema: "https://proj.org/schemas/v0.7/projjson.schema.json", + type: "GeographicCRS", + name: "WGS 84 (CRS84)", + datum_ensemble: { + name: "World Geodetic System 1984 ensemble", + members: [ + { + name: "World Geodetic System 1984 (Transit)", + id: { authority: "EPSG", code: 1166 }, + }, + { + name: "World Geodetic System 1984 (G730)", + id: { authority: "EPSG", code: 1152 }, + }, + { + name: "World Geodetic System 1984 (G873)", + id: { authority: "EPSG", code: 1153 }, + }, + { + name: "World Geodetic System 1984 (G1150)", + id: { authority: "EPSG", code: 1154 }, + }, + { + name: "World Geodetic System 1984 (G1674)", + id: { authority: "EPSG", code: 1155 }, + }, + { + name: "World Geodetic System 1984 (G1762)", + id: { authority: "EPSG", code: 1156 }, + }, + { + name: "World Geodetic System 1984 (G2139)", + id: { authority: "EPSG", code: 1309 }, + }, + ], + ellipsoid: { + name: "WGS 84", + semi_major_axis: 6378137, + inverse_flattening: 298.257223563, + }, + accuracy: "2.0", + id: { authority: "EPSG", code: 6326 }, + }, + coordinate_system: { + subtype: "ellipsoidal", + axis: [ + { + name: "Geodetic longitude", + abbreviation: "Lon", + direction: "east", + unit: "degree", + }, + { + name: "Geodetic latitude", + abbreviation: "Lat", + direction: "north", + unit: "degree", + }, + ], + }, + scope: "Not known.", + area: "World.", + bbox: { + south_latitude: -90, + west_longitude: -180, + north_latitude: 90, + east_longitude: 180, + }, + id: { authority: "OGC", code: "CRS84" }, +}; + +/** + * Extract COG metadata + */ +export async function extractCOGMetadata(tiff: GeoTIFF): Promise { + const image = await tiff.getImage(); + + const width = image.getWidth(); + const height = image.getHeight(); + const tileWidth = image.getTileWidth(); + const tileHeight = image.getTileHeight(); + + const tilesX = Math.ceil(width / tileWidth); + const tilesY = Math.ceil(height / tileHeight); + + const bbox = image.getBoundingBox(); + const geoKeys = image.getGeoKeys(); + const projectionCode: number | null = + geoKeys.ProjectedCSTypeGeoKey || geoKeys.GeographicTypeGeoKey || null; + const projection = projectionCode ? `EPSG:${projectionCode}` : null; + + // Extract geotransform from full-resolution image + // Only the top-level IFD has geo keys, so we'll derive overviews from this + const baseGeotransform = extractGeotransform(image); + + // Overviews **in COG order**, from finest to coarsest (we'll reverse the + // array later) + const overviews: COGOverview[] = []; + const imageCount = await tiff.getImageCount(); + + // Full resolution image (GeoTIFF index 0) + overviews.push({ + geoTiffIndex: 0, + width, + height, + tilesX, + tilesY, + scaleFactor: 1, + geotransform: baseGeotransform, + // TODO: combine these two properties into one + level: imageCount - 1, // Coarsest level number + z: imageCount - 1, + }); + + for (let i = 1; i < imageCount; i++) { + const overview = await tiff.getImage(i); + const overviewWidth = overview.getWidth(); + const overviewHeight = overview.getHeight(); + const overviewTileWidth = overview.getTileWidth(); + const overviewTileHeight = overview.getTileHeight(); + + const scaleFactor = Math.round(width / overviewWidth); + + // Derive geotransform for this overview by scaling pixel size + // [a, b, c, d, e, f] where a and e are pixel dimensions + const overviewGeotransform: [ + number, + number, + number, + number, + number, + number, + ] = [ + baseGeotransform[0] * scaleFactor, // a: scaled pixel width + baseGeotransform[1] * scaleFactor, // b: scaled row rotation + baseGeotransform[2], // c: same x origin + baseGeotransform[3] * scaleFactor, // d: scaled column rotation + baseGeotransform[4] * scaleFactor, // e: scaled pixel height (typically negative) + baseGeotransform[5], // f: same y origin + ]; + + overviews.push({ + geoTiffIndex: i, + width: overviewWidth, + height: overviewHeight, + tilesX: Math.ceil(overviewWidth / overviewTileWidth), + tilesY: Math.ceil(overviewHeight / overviewTileHeight), + scaleFactor, + geotransform: overviewGeotransform, + // TODO: combine these two properties into one + level: imageCount - 1 - i, + z: imageCount - 1 - i, + }); + } + + // Reverse to TileMatrixSet order: coarsest (0) → finest (n) + overviews.reverse(); + + const sourceProjection = await getProjjson(projectionCode); + const projectToWgs84 = proj4(sourceProjection, OGC_84); + const projectTo3857 = proj4(sourceProjection, "EPSG:3857"); + + return { + width, + height, + tileWidth, + tileHeight, + tilesX, + tilesY, + bbox: [bbox[0], bbox[1], bbox[2], bbox[3]], + projection, + projectToWgs84, + projectTo3857, + overviews, + image: tiff, + }; +} + +async function getProjjson(projectionCode: number | null) { + const url = `https://epsg.io/${projectionCode}.json`; + const response = await fetch(url); + if (!response.ok) { + throw new Error(`Failed to fetch projection data from ${url}`); + } + const data = await response.json(); + return data; +} + +const viewport = new WebMercatorViewport({ + height: 500, + width: 845, + latitude: 40.88775942857086, + longitude: -73.20197979318772, + zoom: 11.294596276534985, +}); + +/** + * COGTileset2D with proper frustum culling + */ +export class COGTileset2D extends Tileset2D { + private cogMetadata: COGMetadata; + + constructor(cogMetadata: COGMetadata, opts: Tileset2DProps) { + super(opts); + this.cogMetadata = cogMetadata; + } + + /** + * Get tile indices visible in viewport + * Uses frustum culling similar to OSM implementation + * + * Overviews follow TileMatrixSet ordering: index 0 = coarsest, higher = finer + */ + getTileIndices(opts: { + viewport: Viewport; + maxZoom?: number; + minZoom?: number; + zRange: ZRange | null; + modelMatrix?: Matrix4; + modelMatrixInverse?: Matrix4; + }): COGTileIndex[] { + console.log("Called getTileIndices", opts); + const tileIndices = getTileIndices(this.cogMetadata, opts); + console.log("Visible tile indices:", tileIndices); + + // return [ + // { x: 0, y: 0, z: 0 }, + // { x: 0, y: 0, z: 1 }, + // { x: 1, y: 1, z: 2 }, + // { x: 1, y: 2, z: 3 }, + // { x: 2, y: 1, z: 3 }, + // { x: 2, y: 2, z: 3 }, + // { x: 3, y: 1, z: 3 }, + // ]; // Temporary override for testing + return tileIndices; + } + + getTileId(index: COGTileIndex): string { + return `${index.x}-${index.y}-${index.z}`; + } + + getParentIndex(index: COGTileIndex): COGTileIndex { + if (index.z === 0) { + // Already at coarsest level + return index; + } + + const currentOverview = this.cogMetadata.overviews[index.z]; + const parentOverview = this.cogMetadata.overviews[index.z - 1]; + + const scaleFactor = + currentOverview.scaleFactor / parentOverview.scaleFactor; + + return { + x: Math.floor(index.x / scaleFactor), + y: Math.floor(index.y / scaleFactor), + z: index.z - 1, + }; + } + + getTileZoom(index: COGTileIndex): number { + return index.z; + } + + getTileMetadata(index: COGTileIndex): Record { + const { x, y, z } = index; + const { overviews, tileWidth, tileHeight } = this.cogMetadata; + const overview = overviews[z]; + + // Use geotransform to calculate tile bounds + // geotransform: [a, b, c, d, e, f] where: + // x_geo = a * col + b * row + c + // y_geo = d * col + e * row + f + const [a, b, c, d, e, f] = overview.geotransform; + + // Calculate pixel coordinates for this tile's extent + const pixelMinCol = x * tileWidth; + const pixelMinRow = y * tileHeight; + const pixelMaxCol = (x + 1) * tileWidth; + const pixelMaxRow = (y + 1) * tileHeight; + + // Calculate the four corners of the tile in geographic coordinates + const topLeft = [ + a * pixelMinCol + b * pixelMinRow + c, + d * pixelMinCol + e * pixelMinRow + f, + ]; + const topRight = [ + a * pixelMaxCol + b * pixelMinRow + c, + d * pixelMaxCol + e * pixelMinRow + f, + ]; + const bottomLeft = [ + a * pixelMinCol + b * pixelMaxRow + c, + d * pixelMinCol + e * pixelMaxRow + f, + ]; + const bottomRight = [ + a * pixelMaxCol + b * pixelMaxRow + c, + d * pixelMaxCol + e * pixelMaxRow + f, + ]; + + // Return the projected bounds as four corners + // This preserves rotation/skew information + const projectedBounds = { + topLeft, + topRight, + bottomLeft, + bottomRight, + }; + + // Also compute axis-aligned bounding box for compatibility + const bounds: Bounds = [ + Math.min(topLeft[0], topRight[0], bottomLeft[0], bottomRight[0]), + Math.min(topLeft[1], topRight[1], bottomLeft[1], bottomRight[1]), + Math.max(topLeft[0], topRight[0], bottomLeft[0], bottomRight[0]), + Math.max(topLeft[1], topRight[1], bottomLeft[1], bottomRight[1]), + ]; + + return { + bounds, + projectedBounds, + tileWidth, + tileHeight, + overview, + }; + } +} From 91d403dd9a9d3d55d2aa521bc16241e7f1a3a600 Mon Sep 17 00:00:00 2001 From: Kyle Barron Date: Mon, 15 Dec 2025 20:17:04 -0500 Subject: [PATCH 02/14] Add dependency --- packages/deck.gl-raster/package.json | 1 + pnpm-lock.yaml | 3 +++ 2 files changed, 4 insertions(+) diff --git a/packages/deck.gl-raster/package.json b/packages/deck.gl-raster/package.json index ff032da..81cb8a9 100644 --- a/packages/deck.gl-raster/package.json +++ b/packages/deck.gl-raster/package.json @@ -59,6 +59,7 @@ }, "dependencies": { "@developmentseed/raster-reproject": "workspace:^", + "@math.gl/core": "^4.1.0", "@math.gl/culling": "^4.1.0" }, "volta": { diff --git a/pnpm-lock.yaml b/pnpm-lock.yaml index a7096c9..8210848 100644 --- a/pnpm-lock.yaml +++ b/pnpm-lock.yaml @@ -133,6 +133,9 @@ importers: '@developmentseed/raster-reproject': specifier: workspace:^ version: link:../raster-reproject + '@math.gl/core': + specifier: ^4.1.0 + version: 4.1.0 '@math.gl/culling': specifier: ^4.1.0 version: 4.1.0 From 2432bab95e217fc4dd735da2e3787750c171fe03 Mon Sep 17 00:00:00 2001 From: Kyle Barron Date: Mon, 15 Dec 2025 20:17:28 -0500 Subject: [PATCH 03/14] remove cog-specific code from generic deck.gl-raster implementation --- .../src/raster-tileset/raster-tileset-2d.ts | 197 +----------------- 1 file changed, 4 insertions(+), 193 deletions(-) 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 e3e8ccf..ad936e3 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 @@ -7,203 +7,14 @@ import { Viewport, WebMercatorViewport } from "@deck.gl/core"; import { _Tileset2D as Tileset2D } from "@deck.gl/geo-layers"; -import type { Tileset2DProps } from "@deck.gl/geo-layers/dist/tileset-2d"; -import type { ZRange } from "@deck.gl/geo-layers/dist/tileset-2d/types"; import { Matrix4 } from "@math.gl/core"; -import { GeoTIFF, GeoTIFFImage } from "geotiff"; -import proj4 from "proj4"; import { getTileIndices } from "./raster-tile-traversal"; -import type { COGMetadata, COGTileIndex, COGOverview, Bounds } from "./types"; +import type { TileIndex, Bounds, TileMatrixSet, ZRange } from "./types"; -const OGC_84 = { - $schema: "https://proj.org/schemas/v0.7/projjson.schema.json", - type: "GeographicCRS", - name: "WGS 84 (CRS84)", - datum_ensemble: { - name: "World Geodetic System 1984 ensemble", - members: [ - { - name: "World Geodetic System 1984 (Transit)", - id: { authority: "EPSG", code: 1166 }, - }, - { - name: "World Geodetic System 1984 (G730)", - id: { authority: "EPSG", code: 1152 }, - }, - { - name: "World Geodetic System 1984 (G873)", - id: { authority: "EPSG", code: 1153 }, - }, - { - name: "World Geodetic System 1984 (G1150)", - id: { authority: "EPSG", code: 1154 }, - }, - { - name: "World Geodetic System 1984 (G1674)", - id: { authority: "EPSG", code: 1155 }, - }, - { - name: "World Geodetic System 1984 (G1762)", - id: { authority: "EPSG", code: 1156 }, - }, - { - name: "World Geodetic System 1984 (G2139)", - id: { authority: "EPSG", code: 1309 }, - }, - ], - ellipsoid: { - name: "WGS 84", - semi_major_axis: 6378137, - inverse_flattening: 298.257223563, - }, - accuracy: "2.0", - id: { authority: "EPSG", code: 6326 }, - }, - coordinate_system: { - subtype: "ellipsoidal", - axis: [ - { - name: "Geodetic longitude", - abbreviation: "Lon", - direction: "east", - unit: "degree", - }, - { - name: "Geodetic latitude", - abbreviation: "Lat", - direction: "north", - unit: "degree", - }, - ], - }, - scope: "Not known.", - area: "World.", - bbox: { - south_latitude: -90, - west_longitude: -180, - north_latitude: 90, - east_longitude: 180, - }, - id: { authority: "OGC", code: "CRS84" }, -}; - -/** - * Extract COG metadata - */ -export async function extractCOGMetadata(tiff: GeoTIFF): Promise { - const image = await tiff.getImage(); - - const width = image.getWidth(); - const height = image.getHeight(); - const tileWidth = image.getTileWidth(); - const tileHeight = image.getTileHeight(); - - const tilesX = Math.ceil(width / tileWidth); - const tilesY = Math.ceil(height / tileHeight); - - const bbox = image.getBoundingBox(); - const geoKeys = image.getGeoKeys(); - const projectionCode: number | null = - geoKeys.ProjectedCSTypeGeoKey || geoKeys.GeographicTypeGeoKey || null; - const projection = projectionCode ? `EPSG:${projectionCode}` : null; - - // Extract geotransform from full-resolution image - // Only the top-level IFD has geo keys, so we'll derive overviews from this - const baseGeotransform = extractGeotransform(image); - - // Overviews **in COG order**, from finest to coarsest (we'll reverse the - // array later) - const overviews: COGOverview[] = []; - const imageCount = await tiff.getImageCount(); - - // Full resolution image (GeoTIFF index 0) - overviews.push({ - geoTiffIndex: 0, - width, - height, - tilesX, - tilesY, - scaleFactor: 1, - geotransform: baseGeotransform, - // TODO: combine these two properties into one - level: imageCount - 1, // Coarsest level number - z: imageCount - 1, - }); - - for (let i = 1; i < imageCount; i++) { - const overview = await tiff.getImage(i); - const overviewWidth = overview.getWidth(); - const overviewHeight = overview.getHeight(); - const overviewTileWidth = overview.getTileWidth(); - const overviewTileHeight = overview.getTileHeight(); - - const scaleFactor = Math.round(width / overviewWidth); - - // Derive geotransform for this overview by scaling pixel size - // [a, b, c, d, e, f] where a and e are pixel dimensions - const overviewGeotransform: [ - number, - number, - number, - number, - number, - number, - ] = [ - baseGeotransform[0] * scaleFactor, // a: scaled pixel width - baseGeotransform[1] * scaleFactor, // b: scaled row rotation - baseGeotransform[2], // c: same x origin - baseGeotransform[3] * scaleFactor, // d: scaled column rotation - baseGeotransform[4] * scaleFactor, // e: scaled pixel height (typically negative) - baseGeotransform[5], // f: same y origin - ]; - - overviews.push({ - geoTiffIndex: i, - width: overviewWidth, - height: overviewHeight, - tilesX: Math.ceil(overviewWidth / overviewTileWidth), - tilesY: Math.ceil(overviewHeight / overviewTileHeight), - scaleFactor, - geotransform: overviewGeotransform, - // TODO: combine these two properties into one - level: imageCount - 1 - i, - z: imageCount - 1 - i, - }); - } - - // Reverse to TileMatrixSet order: coarsest (0) → finest (n) - overviews.reverse(); - - const sourceProjection = await getProjjson(projectionCode); - const projectToWgs84 = proj4(sourceProjection, OGC_84); - const projectTo3857 = proj4(sourceProjection, "EPSG:3857"); - - return { - width, - height, - tileWidth, - tileHeight, - tilesX, - tilesY, - bbox: [bbox[0], bbox[1], bbox[2], bbox[3]], - projection, - projectToWgs84, - projectTo3857, - overviews, - image: tiff, - }; -} - -async function getProjjson(projectionCode: number | null) { - const url = `https://epsg.io/${projectionCode}.json`; - const response = await fetch(url); - if (!response.ok) { - throw new Error(`Failed to fetch projection data from ${url}`); - } - const data = await response.json(); - return data; -} +// Include correct Tileset2DProps type when exported from deck.gl +// https://github.com/visgl/deck.gl/pull/9917 +type Tileset2DProps = any; const viewport = new WebMercatorViewport({ height: 500, From db41bde9f0e4f8140774da63331d9c3c24748f45 Mon Sep 17 00:00:00 2001 From: Kyle Barron Date: Mon, 15 Dec 2025 20:18:28 -0500 Subject: [PATCH 04/14] type renames --- .../src/raster-tileset/raster-tileset-2d.ts | 42 ++++++++----------- 1 file changed, 17 insertions(+), 25 deletions(-) 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 ad936e3..0ad672e 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 @@ -1,11 +1,11 @@ /** - * COGTileset2D - Improved Implementation with Frustum Culling + * Rasterileset2D - Improved Implementation with Frustum Culling * * This version properly implements frustum culling and bounding volume calculations * following the pattern from deck.gl's OSM tile indexing. */ -import { Viewport, WebMercatorViewport } from "@deck.gl/core"; +import { Viewport } from "@deck.gl/core"; import { _Tileset2D as Tileset2D } from "@deck.gl/geo-layers"; import { Matrix4 } from "@math.gl/core"; @@ -16,23 +16,15 @@ import type { TileIndex, Bounds, TileMatrixSet, ZRange } from "./types"; // https://github.com/visgl/deck.gl/pull/9917 type Tileset2DProps = any; -const viewport = new WebMercatorViewport({ - height: 500, - width: 845, - latitude: 40.88775942857086, - longitude: -73.20197979318772, - zoom: 11.294596276534985, -}); - /** - * COGTileset2D with proper frustum culling + * Rasterileset2D with proper frustum culling */ -export class COGTileset2D extends Tileset2D { - private cogMetadata: COGMetadata; +export class Rasterileset2D extends Tileset2D { + private metadata: TileMatrixSet; - constructor(cogMetadata: COGMetadata, opts: Tileset2DProps) { + constructor(metadata: TileMatrixSet, opts: Tileset2DProps) { super(opts); - this.cogMetadata = cogMetadata; + this.metadata = metadata; } /** @@ -41,16 +33,16 @@ export class COGTileset2D extends Tileset2D { * * Overviews follow TileMatrixSet ordering: index 0 = coarsest, higher = finer */ - getTileIndices(opts: { + override getTileIndices(opts: { viewport: Viewport; maxZoom?: number; minZoom?: number; zRange: ZRange | null; modelMatrix?: Matrix4; modelMatrixInverse?: Matrix4; - }): COGTileIndex[] { + }): TileIndex[] { console.log("Called getTileIndices", opts); - const tileIndices = getTileIndices(this.cogMetadata, opts); + const tileIndices = getTileIndices(this.metadata, opts); console.log("Visible tile indices:", tileIndices); // return [ @@ -65,18 +57,18 @@ export class COGTileset2D extends Tileset2D { return tileIndices; } - getTileId(index: COGTileIndex): string { + override getTileId(index: TileIndex): string { return `${index.x}-${index.y}-${index.z}`; } - getParentIndex(index: COGTileIndex): COGTileIndex { + override getParentIndex(index: TileIndex): TileIndex { if (index.z === 0) { // Already at coarsest level return index; } - const currentOverview = this.cogMetadata.overviews[index.z]; - const parentOverview = this.cogMetadata.overviews[index.z - 1]; + const currentOverview = this.metadata.tileMatrices[index.z]; + const parentOverview = this.metadata.tileMatrices[index.z - 1]; const scaleFactor = currentOverview.scaleFactor / parentOverview.scaleFactor; @@ -88,13 +80,13 @@ export class COGTileset2D extends Tileset2D { }; } - getTileZoom(index: COGTileIndex): number { + override getTileZoom(index: TileIndex): number { return index.z; } - getTileMetadata(index: COGTileIndex): Record { + override getTileMetadata(index: TileIndex): Record { const { x, y, z } = index; - const { overviews, tileWidth, tileHeight } = this.cogMetadata; + const { overviews, tileWidth, tileHeight } = this.metadata; const overview = overviews[z]; // Use geotransform to calculate tile bounds From 9a8d2df1861d779ae44bcede673513b725d7ab8e Mon Sep 17 00:00:00 2001 From: Kyle Barron Date: Mon, 15 Dec 2025 21:26:08 -0500 Subject: [PATCH 05/14] edit --- packages/deck.gl-raster/src/raster-tileset/index.ts | 2 ++ .../deck.gl-raster/src/raster-tileset/raster-tileset-2d.ts | 6 +++--- 2 files changed, 5 insertions(+), 3 deletions(-) create mode 100644 packages/deck.gl-raster/src/raster-tileset/index.ts diff --git a/packages/deck.gl-raster/src/raster-tileset/index.ts b/packages/deck.gl-raster/src/raster-tileset/index.ts new file mode 100644 index 0000000..d19311b --- /dev/null +++ b/packages/deck.gl-raster/src/raster-tileset/index.ts @@ -0,0 +1,2 @@ +export { RasterTileset2D } from "./raster-tileset-2d.js"; +export type { TileMatrix, TileMatrixSet } from "./types.js"; 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 0ad672e..12bb20e 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 @@ -1,5 +1,5 @@ /** - * Rasterileset2D - Improved Implementation with Frustum Culling + * RasterTileset2D - Improved Implementation with Frustum Culling * * This version properly implements frustum culling and bounding volume calculations * following the pattern from deck.gl's OSM tile indexing. @@ -17,9 +17,9 @@ import type { TileIndex, Bounds, TileMatrixSet, ZRange } from "./types"; type Tileset2DProps = any; /** - * Rasterileset2D with proper frustum culling + * RasterTileset2D with proper frustum culling */ -export class Rasterileset2D extends Tileset2D { +export class RasterTileset2D extends Tileset2D { private metadata: TileMatrixSet; constructor(metadata: TileMatrixSet, opts: Tileset2DProps) { From 9b1a667af61b6c467f1d94788a10b96ca044d31a Mon Sep 17 00:00:00 2001 From: Kyle Barron Date: Mon, 15 Dec 2025 21:33:15 -0500 Subject: [PATCH 06/14] cleanup --- examples/cog-basic/package.json | 1 + packages/deck.gl-raster/src/index.ts | 3 +- .../src/raster-tileset/raster-tileset-2d.ts | 28 +++++++++---------- pnpm-lock.yaml | 3 ++ 4 files changed, 20 insertions(+), 15 deletions(-) diff --git a/examples/cog-basic/package.json b/examples/cog-basic/package.json index ebcf94b..9ecf444 100644 --- a/examples/cog-basic/package.json +++ b/examples/cog-basic/package.json @@ -14,6 +14,7 @@ "@deck.gl/mapbox": "^9.2.5", "@deck.gl/mesh-layers": "^9.2.5", "@developmentseed/deck.gl-cog": "workspace:^", + "@developmentseed/deck.gl-raster": "workspace:^", "geotiff": "^2.1.3", "maplibre-gl": "^5.0.0", "proj4": "^2.20.2", diff --git a/packages/deck.gl-raster/src/index.ts b/packages/deck.gl-raster/src/index.ts index 47f8466..46a1518 100644 --- a/packages/deck.gl-raster/src/index.ts +++ b/packages/deck.gl-raster/src/index.ts @@ -1,3 +1,4 @@ export { RasterLayer } from "./raster-layer.js"; export type { RasterLayerProps } from "./raster-layer.js"; -export type { TileMatrixSet, TileMatrix } from "./raster-tileset/types.js"; +export { RasterTileset2D } from "./raster-tileset/index.js"; +export type { TileMatrix, TileMatrixSet } from "./raster-tileset/types.js"; 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 12bb20e..1fa72b8 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 @@ -67,15 +67,14 @@ export class RasterTileset2D extends Tileset2D { return index; } - const currentOverview = this.metadata.tileMatrices[index.z]; - const parentOverview = this.metadata.tileMatrices[index.z - 1]; + const currentOverview = this.metadata.tileMatrices[index.z]!; + const parentOverview = this.metadata.tileMatrices[index.z - 1]!; - const scaleFactor = - currentOverview.scaleFactor / parentOverview.scaleFactor; + const decimation = currentOverview.cellSize / parentOverview.cellSize; return { - x: Math.floor(index.x / scaleFactor), - y: Math.floor(index.y / scaleFactor), + x: Math.floor(index.x / decimation), + y: Math.floor(index.y / decimation), z: index.z - 1, }; } @@ -86,14 +85,15 @@ export class RasterTileset2D extends Tileset2D { override getTileMetadata(index: TileIndex): Record { const { x, y, z } = index; - const { overviews, tileWidth, tileHeight } = this.metadata; - const overview = overviews[z]; + const { tileMatrices } = this.metadata; + const tileMatrix = tileMatrices[z]!; + const { geotransform, tileHeight, tileWidth } = tileMatrix; // Use geotransform to calculate tile bounds // geotransform: [a, b, c, d, e, f] where: // x_geo = a * col + b * row + c // y_geo = d * col + e * row + f - const [a, b, c, d, e, f] = overview.geotransform; + const [a, b, c, d, e, f] = geotransform; // Calculate pixel coordinates for this tile's extent const pixelMinCol = x * tileWidth; @@ -102,19 +102,19 @@ export class RasterTileset2D extends Tileset2D { const pixelMaxRow = (y + 1) * tileHeight; // Calculate the four corners of the tile in geographic coordinates - const topLeft = [ + const topLeft: [number, number] = [ a * pixelMinCol + b * pixelMinRow + c, d * pixelMinCol + e * pixelMinRow + f, ]; - const topRight = [ + const topRight: [number, number] = [ a * pixelMaxCol + b * pixelMinRow + c, d * pixelMaxCol + e * pixelMinRow + f, ]; - const bottomLeft = [ + const bottomLeft: [number, number] = [ a * pixelMinCol + b * pixelMaxRow + c, d * pixelMinCol + e * pixelMaxRow + f, ]; - const bottomRight = [ + const bottomRight: [number, number] = [ a * pixelMaxCol + b * pixelMaxRow + c, d * pixelMaxCol + e * pixelMaxRow + f, ]; @@ -141,7 +141,7 @@ export class RasterTileset2D extends Tileset2D { projectedBounds, tileWidth, tileHeight, - overview, + tileMatrix, }; } } diff --git a/pnpm-lock.yaml b/pnpm-lock.yaml index 8210848..180517a 100644 --- a/pnpm-lock.yaml +++ b/pnpm-lock.yaml @@ -35,6 +35,9 @@ importers: '@developmentseed/deck.gl-cog': specifier: workspace:^ version: link:../../packages/deck.gl-cog + '@developmentseed/deck.gl-raster': + specifier: workspace:^ + version: link:../../packages/deck.gl-raster geotiff: specifier: ^2.1.3 version: 2.1.3 From 2fc72aad02097dc9b954daf4d7b8438d08f06a85 Mon Sep 17 00:00:00 2001 From: Kyle Barron Date: Tue, 16 Dec 2025 11:34:13 -0500 Subject: [PATCH 07/14] add dep --- examples/cog-basic/package.json | 1 + pnpm-lock.yaml | 3 +++ 2 files changed, 4 insertions(+) diff --git a/examples/cog-basic/package.json b/examples/cog-basic/package.json index 9ecf444..48fe143 100644 --- a/examples/cog-basic/package.json +++ b/examples/cog-basic/package.json @@ -11,6 +11,7 @@ "dependencies": { "@deck.gl/core": "^9.2.5", "@deck.gl/layers": "^9.2.5", + "@deck.gl/geo-layers": "^9.2.5", "@deck.gl/mapbox": "^9.2.5", "@deck.gl/mesh-layers": "^9.2.5", "@developmentseed/deck.gl-cog": "workspace:^", diff --git a/pnpm-lock.yaml b/pnpm-lock.yaml index 180517a..dca5a65 100644 --- a/pnpm-lock.yaml +++ b/pnpm-lock.yaml @@ -23,6 +23,9 @@ importers: '@deck.gl/core': specifier: ^9.2.5 version: 9.2.5 + '@deck.gl/geo-layers': + specifier: ^9.2.5 + version: 9.2.5(@deck.gl/core@9.2.5)(@deck.gl/extensions@9.2.5(@deck.gl/core@9.2.5)(@luma.gl/core@9.2.4)(@luma.gl/engine@9.2.4(@luma.gl/core@9.2.4)(@luma.gl/shadertools@9.2.4(@luma.gl/core@9.2.4))))(@deck.gl/layers@9.2.5(@deck.gl/core@9.2.5)(@loaders.gl/core@4.3.4)(@luma.gl/core@9.2.4)(@luma.gl/engine@9.2.4(@luma.gl/core@9.2.4)(@luma.gl/shadertools@9.2.4(@luma.gl/core@9.2.4))))(@deck.gl/mesh-layers@9.2.5(@deck.gl/core@9.2.5)(@loaders.gl/core@4.3.4)(@luma.gl/core@9.2.4)(@luma.gl/engine@9.2.4(@luma.gl/core@9.2.4)(@luma.gl/shadertools@9.2.4(@luma.gl/core@9.2.4)))(@luma.gl/gltf@9.2.4(@luma.gl/constants@9.2.4)(@luma.gl/core@9.2.4)(@luma.gl/engine@9.2.4(@luma.gl/core@9.2.4)(@luma.gl/shadertools@9.2.4(@luma.gl/core@9.2.4)))(@luma.gl/shadertools@9.2.4(@luma.gl/core@9.2.4)))(@luma.gl/shadertools@9.2.4(@luma.gl/core@9.2.4)))(@loaders.gl/core@4.3.4)(@luma.gl/constants@9.2.4)(@luma.gl/core@9.2.4)(@luma.gl/engine@9.2.4(@luma.gl/core@9.2.4)(@luma.gl/shadertools@9.2.4(@luma.gl/core@9.2.4))) '@deck.gl/layers': specifier: ^9.2.5 version: 9.2.5(@deck.gl/core@9.2.5)(@loaders.gl/core@4.3.4)(@luma.gl/core@9.2.4)(@luma.gl/engine@9.2.4(@luma.gl/core@9.2.4)(@luma.gl/shadertools@9.2.4(@luma.gl/core@9.2.4))) From 7e96986e2f0ea4d0951670aaa57510850bd4f82a Mon Sep 17 00:00:00 2001 From: Kyle Barron Date: Tue, 16 Dec 2025 11:34:44 -0500 Subject: [PATCH 08/14] wip debug tile layer in app --- examples/cog-basic/src/App.tsx | 106 +++++++++++++++++++++++++++++---- 1 file changed, 94 insertions(+), 12 deletions(-) diff --git a/examples/cog-basic/src/App.tsx b/examples/cog-basic/src/App.tsx index 06ad9fd..99b4254 100644 --- a/examples/cog-basic/src/App.tsx +++ b/examples/cog-basic/src/App.tsx @@ -1,13 +1,22 @@ import { useEffect, useState, useRef } from "react"; import { Map, useControl, type MapRef } from "react-map-gl/maplibre"; +import type { Tileset2DProps } from "@deck.gl/geo-layers/dist/tileset-2d"; import { MapboxOverlay } from "@deck.gl/mapbox"; +import { PathLayer } from "@deck.gl/layers"; import type { DeckProps } from "@deck.gl/core"; +import { TileLayer, TileLayerProps } from "@deck.gl/geo-layers"; import { fromUrl } from "geotiff"; import type { GeoTIFF } from "geotiff"; -import { COGLayer } from "@developmentseed/deck.gl-cog"; +import { COGLayer, parseCOGTileMatrixSet } from "@developmentseed/deck.gl-cog"; +import { + RasterTileset2D, + TileMatrixSet, +} from "@developmentseed/deck.gl-raster"; import proj4 from "proj4"; import "maplibre-gl/dist/maplibre-gl.css"; +window.proj4 = proj4; + function DeckGLOverlay(props: DeckProps) { const overlay = useControl(() => new MapboxOverlay(props)); overlay.setProps(props); @@ -74,6 +83,7 @@ const COG_URL = export default function App() { const mapRef = useRef(null); const [geotiff, setGeotiff] = useState(null); + const [cogMetadata, setCogMetadata] = useState(null); const [error, setError] = useState(null); const [loading, setLoading] = useState(true); const [debug, setDebug] = useState(false); @@ -88,10 +98,17 @@ export default function App() { setError(null); const tiff = await fromUrl(COG_URL); + window.tiff = tiff; if (mounted) { setGeotiff(tiff); + const m = await parseCOGTileMatrixSet(tiff); + console.log("COG TileMatrixSet:", m); + window.m = m; + setCogMetadata(m); + // window.cogMetadata = cogMetadata;d + // Calculate bounds and fit to them const bounds = await getCogBounds(tiff); if (mapRef.current) { @@ -120,17 +137,22 @@ export default function App() { }; }, []); - const layers = geotiff - ? [ - new COGLayer({ - id: "cog-layer", - geotiff, - maxError: 0.125, - debug, - debugOpacity, - }), - ] - : []; + const layers = + geotiff && cogMetadata + ? [ + // new COGLayer({ + // id: "cog-layer", + // geotiff, + // maxError: 0.125, + // debug, + // debugOpacity, + // }), + createTileLayer(cogMetadata, { + id: "raster-tile-layer", + data: geotiff, + }), + ] + : []; return (
@@ -286,3 +308,63 @@ export default function App() {
); } + +function createTileLayer(metadata: TileMatrixSet, props: TileLayerProps) { + // Create a factory class that wraps COGTileset2D with the metadata + class RasterTilesetWrapper extends RasterTileset2D { + constructor(opts: Tileset2DProps) { + super(metadata, opts); + } + } + + return new TileLayer({ + ...props, + TilesetClass: RasterTilesetWrapper, + renderSubLayers: (props) => { + const { tile } = props; + console.log("Rendering tile:", tile); + + // Get projected bounds from tile data + // getTileMetadata returns data that includes projectedBounds + // eslint-disable-next-line @typescript-eslint/no-explicit-any + const projectedBounds = (tile as any)?.projectedBounds; + + if (!projectedBounds || !metadata) { + return []; + } + + // Project bounds from image CRS to WGS84 + const { topLeft, topRight, bottomLeft, bottomRight } = projectedBounds; + + const topLeftWgs84 = metadata.projectToWgs84(topLeft); + const topRightWgs84 = metadata.projectToWgs84(topRight); + const bottomRightWgs84 = metadata.projectToWgs84(bottomRight); + const bottomLeftWgs84 = metadata.projectToWgs84(bottomLeft); + + // Create a closed path around the tile bounds + const path = [ + topLeftWgs84, + topRightWgs84, + bottomRightWgs84, + bottomLeftWgs84, + topLeftWgs84, // Close the path + ]; + + console.log("Tile bounds path (WGS84):", path); + + return [ + new PathLayer({ + id: `${tile.id}-bounds`, + data: [{ path }], + getPath: (d) => d.path, + getColor: [255, 0, 0, 255], // Red + getWidth: 2, + widthUnits: "pixels", + pickable: false, + }), + ]; + + // return null; + }, + }); +} From ec032d2b05b588ebcd8b86965d08aa3a01d7d981 Mon Sep 17 00:00:00 2001 From: Kyle Barron Date: Tue, 16 Dec 2025 12:26:46 -0500 Subject: [PATCH 09/14] progress on close inspection of traversal --- .../raster-tileset/raster-tile-traversal.ts | 66 +++++++++++-------- 1 file changed, 38 insertions(+), 28 deletions(-) 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 c6f17fb..8ecc171 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 @@ -1,8 +1,9 @@ /** - * This file implements tile traversal for generic 2D tilesets defined by COG - * tile layouts. + * This file implements tile traversal for generic 2D tilesets defined by + * TileMatrixSet tile layouts. * * The main algorithm works as follows: + * * 1. Start at the root tile(s) (z=0, covers the entire image, but not * necessarily the whole world) * 2. Test if each tile is visible using viewport frustum culling @@ -22,14 +23,14 @@ import { } from "@deck.gl/core"; import { CullingVolume, - Plane, makeOrientedBoundingBoxFromPoints, + Plane, } from "@math.gl/culling"; import type { - TileMatrixSet, - TileMatrix, TileIndex, + TileMatrix, + TileMatrixSet, ZRange, } from "../raster-tileset/types.js"; @@ -38,10 +39,13 @@ import type { * * The world always spans [0, 512] in both X and Y in Web Mercator common space. * + * At zoom level 0, there is 1 tile that represents the whole world, so that tile is 512x512 units. + * At zoom level z, there are 2^z tiles along each axis, so each tile is (512 / 2^z) units. + * * The origin (0,0) is at the top-left corner, and (512,512) is at the * bottom-right. */ -const WORLD_SIZE = 512; +const TILE_SIZE = 512; // Reference points used to sample tile boundaries for bounding volume // calculation. @@ -78,23 +82,24 @@ const REF_POINTS_9 = REF_POINTS_5.concat([ ]); /** - * Raster Tile Node - similar to OSMNode but for a generic raster tileset's - * tile structure. + * Raster Tile Node - represents a single tile in the TileMatrixSet structure * - * Represents a single tile in the COG internal tiling pyramid. + * Akin to the upstream OSMNode class. * - * COG tile nodes use the following coordinate system: + * This node class uses the following coordinate system: * - * - x: tile column (0 to TileMatrix.tilesX, left to right) - * - y: tile row (0 to TileMatrix.tilesY, top to bottom) - * - z: overview level. This uses TileMatrixSet ordering where: 0 = coarsest, higher = finer + * - x: tile column (0 to TileMatrix.matrixWidth, left to right) + * - y: tile row (0 to TileMatrix.matrixHeight, top to bottom) + * - z: overview level. This assumes ordering where: 0 = coarsest, higher = finer */ export class RasterTileNode { /** Index across a row */ x: number; + /** Index down a column */ y: number; - /** TileMatrixSet-style zoom index (higher = finer detail) */ + + /** Zoom index assumed to be (higher = finer detail) */ z: number; private metadata: TileMatrixSet; @@ -125,40 +130,45 @@ export class RasterTileNode { } /** Get overview info for this tile's z level */ - get overview(): TileMatrix { + get tileMatrix(): TileMatrix { return this.metadata.tileMatrices[this.z]!; } /** Get the children of this node. */ - get children(): RasterTileNode[] { + get children(): RasterTileNode[] | null { if (!this._children) { const maxZ = this.metadata.tileMatrices.length - 1; if (this.z >= maxZ) { // Already at finest resolution, no children - return []; + return null; } // In TileMatrixSet ordering: refine to z + 1 (finer detail) const childZ = this.z + 1; - const parentOverview = this.overview; - const childOverview = this.metadata.tileMatrices[childZ]!; - - // Calculate scale factor between levels - const scaleFactor = parentOverview.cellSize / childOverview.cellSize; + const parentMatrix = this.tileMatrix; + const childMatrix = this.metadata.tileMatrices[childZ]!; + + // Calculate decimation between levels + // Note: here we assume that the decimation is an integer. + // For non-integer decimation, the tile origin wouldn't necessarily be in + // the same place as its children. + const decimation = Math.round( + parentMatrix.cellSize / childMatrix.cellSize, + ); // Generate child tiles this._children = []; - for (let dy = 0; dy < scaleFactor; dy++) { - for (let dx = 0; dx < scaleFactor; dx++) { - const childX = this.x * scaleFactor + dx; - const childY = this.y * scaleFactor + dy; + for (let dy = 0; dy < decimation; dy++) { + for (let dx = 0; dx < decimation; dx++) { + const childX = this.x * decimation + dx; + const childY = this.y * decimation + dy; // Only create child if it's within bounds // Some tiles on the edges might not need to be created at higher // resolutions (higher map zoom level) if ( - childX < childOverview.tileWidth && - childY < childOverview.tileHeight + childX < childMatrix.matrixWidth && + childY < childMatrix.matrixHeight ) { this._children.push( new RasterTileNode(childX, childY, childZ, this.metadata), From dee6a31918bb77a2922e85e661e459dd3fb19d3b Mon Sep 17 00:00:00 2001 From: Kyle Barron Date: Tue, 16 Dec 2025 12:52:37 -0500 Subject: [PATCH 10/14] Compute projected tile bounds --- .../raster-tileset/raster-tile-traversal.ts | 88 ++++++++++++++----- 1 file changed, 65 insertions(+), 23 deletions(-) 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 8ecc171..c68ebfd 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 @@ -328,27 +328,27 @@ export class RasterTileNode { /** * Calculate the 3D bounding volume for this tile in deck.gl's common * coordinate space for frustum culling. - * */ getBoundingVolume( zRange: ZRange, project: ((xyz: number[]) => number[]) | null, ) { - const overview = this.overview; - const { tileMatrices } = this.metadata; - const { tileWidth, tileHeight } = tileMatrices[tileMatrices.length - 1]!; - - // Use geotransform to calculate tile bounds - // geotransform: [a, b, c, d, e, f] where: - // x_geo = a * col + b * row + c - // y_geo = d * col + e * row + f - const [a, b, c, d, e, f] = overview.geotransform; - - // Calculate pixel coordinates for this tile's extent - const pixelMinCol = this.x * tileWidth; - const pixelMinRow = this.y * tileHeight; - const pixelMaxCol = (this.x + 1) * tileWidth; - const pixelMaxRow = (this.y + 1) * tileHeight; + const tileMatrix = this.tileMatrix; + const { tileWidth, tileHeight, geotransform } = tileMatrix; + + const projectedBounds = computeProjectedTileBounds({ + x: this.x, + y: this.y, + transform: geotransform, + tileWidth, + tileHeight, + }); + + // I have: + // - a tile index x, y, z + // - tileMatrix is at the current z + // + // I need first the projected bounds // Sample reference points across the tile surface const refPoints = REF_POINTS_9; @@ -451,15 +451,57 @@ export class RasterTileNode { return obb; } +} - /** - * Convert COG coordinates to lng/lat - * This is a placeholder - needs proper projection library (proj4js) - */ - private cogCoordsToLngLat([x, y]: [number, number]): number[] { - const [lng, lat] = this.metadata.projectToWgs84([x, y]); - return [lng, lat, 0]; +/** + * Compute the projected tile bounds in the tile matrix's CRS. + * + * Because it's a linear transformation from the tile index to projected bounds, + * we don't need to sample this for each of the reference points. We only need + * the corners. + * + * @return The bounding box as [minX, minY, maxX, maxY] in projected CRS. + */ +function computeProjectedTileBounds({ + x, + y, + transform, + tileWidth, + tileHeight, +}: { + x: number; + y: number; + transform: [number, number, number, number, number, number]; + tileWidth: number; + tileHeight: number; +}): [number, number, number, number] { + // geotransform: [a, b, c, d, e, f] where: + // x_geo = a * col + b * row + c + // y_geo = d * col + e * row + f + const [a, b, c, d, e, f] = transform; + + // Currently only support non-rotated/non-skewed transforms + if (b !== 0 || d !== 0) { + throw new Error( + `Rotated/skewed geotransforms not yet supported (b=${b}, d=${d}). ` + + `Only north-up, non-rotated rasters are currently supported.`, + ); } + + // Calculate pixel coordinates for this tile's extent + const pixelMinCol = x * tileWidth; + const pixelMinRow = y * tileHeight; + const pixelMaxCol = (x + 1) * tileWidth; + const pixelMaxRow = (y + 1) * tileHeight; + + // Convert pixel coordinates to geographic coordinates using geotransform + const minX = a * pixelMinCol + b * pixelMinRow + c; + const minY = d * pixelMinCol + e * pixelMinRow + f; + + const maxX = a * pixelMaxCol + b * pixelMaxRow + c; + const maxY = d * pixelMaxCol + e * pixelMaxRow + f; + + return [minX, minY, maxX, maxY]; } /** From 05491f0b390e0bb6e592f04b3a42739a83a1951f Mon Sep 17 00:00:00 2001 From: Kyle Barron Date: Tue, 16 Dec 2025 12:53:01 -0500 Subject: [PATCH 11/14] apply rename and edit commetn --- .../raster-tileset/raster-tile-traversal.ts | 25 +++++++++++++------ 1 file changed, 18 insertions(+), 7 deletions(-) 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 c68ebfd..2fd5d4e 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 @@ -181,13 +181,24 @@ export class RasterTileNode { } /** - * Update tile visibility using frustum culling - * This follows the pattern from OSMNode + * Recursively traverse the tile quadtree to determine if this tile (or its descendants) + * should be rendered. + * + * The algorithm performs: + * 1. Visibility culling - reject tiles outside the view frustum + * 2. Bounds checking - reject tiles outside the specified geographic bounds + * 3. LOD selection - choose appropriate zoom level based on distance from camera + * 4. Recursive subdivision - if LOD is insufficient, test child tiles + * + * @returns true if this tile or any descendant is visible, false otherwise */ update(params: { viewport: Viewport; + // Projection: [lng,lat,z] -> common space. Null for Web Mercator. project: ((xyz: number[]) => number[]) | null; + // Camera frustum for visibility testing cullingVolume: CullingVolume; + // [min, max] elevation in common space elevationBounds: ZRange; /** Minimum (coarsest) COG overview level */ minZ: number; @@ -414,12 +425,12 @@ export class RasterTileNode { for (const [mercX, mercY] of refPointPositionsProjected) { // X: offset from [-20M, 20M] to [0, 40M], then normalize to [0, 512] const worldX = - ((mercX + WEB_MERCATOR_MAX) / (2 * WEB_MERCATOR_MAX)) * WORLD_SIZE; + ((mercX + WEB_MERCATOR_MAX) / (2 * WEB_MERCATOR_MAX)) * TILE_SIZE; // Y: same transformation WITHOUT flip // Testing hypothesis: Y-flip might be incorrect since geotransform already handles orientation const worldY = - ((mercY + WEB_MERCATOR_MAX) / (2 * WEB_MERCATOR_MAX)) * WORLD_SIZE; + ((mercY + WEB_MERCATOR_MAX) / (2 * WEB_MERCATOR_MAX)) * TILE_SIZE; console.log( `WebMerc [${mercX.toFixed(2)}, ${mercY.toFixed(2)}] -> World [${worldX.toFixed(4)}, ${worldY.toFixed(4)}]`, @@ -433,10 +444,10 @@ export class RasterTileNode { if (zRange[0] !== zRange[1]) { for (const [mercX, mercY] of refPointPositionsProjected) { const worldX = - ((mercX + WEB_MERCATOR_MAX) / (2 * WEB_MERCATOR_MAX)) * WORLD_SIZE; + ((mercX + WEB_MERCATOR_MAX) / (2 * WEB_MERCATOR_MAX)) * TILE_SIZE; const worldY = - WORLD_SIZE - - ((mercY + WEB_MERCATOR_MAX) / (2 * WEB_MERCATOR_MAX)) * WORLD_SIZE; + TILE_SIZE - + ((mercY + WEB_MERCATOR_MAX) / (2 * WEB_MERCATOR_MAX)) * TILE_SIZE; refPointPositionsWorld.push([worldX, worldY, zRange[1]]); } From 12e350ec9302f2a7220471a5c8f1708e49f172d9 Mon Sep 17 00:00:00 2001 From: Kyle Barron Date: Tue, 16 Dec 2025 14:12:27 -0500 Subject: [PATCH 12/14] First steps of generic bounding volume computation --- .../raster-tileset/raster-tile-traversal.ts | 135 ++++++++++++------ 1 file changed, 89 insertions(+), 46 deletions(-) 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 2fd5d4e..9f9a347 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 @@ -24,6 +24,7 @@ import { import { CullingVolume, makeOrientedBoundingBoxFromPoints, + OrientedBoundingBox, Plane, } from "@math.gl/culling"; @@ -339,61 +340,33 @@ export class RasterTileNode { /** * Calculate the 3D bounding volume for this tile in deck.gl's common * coordinate space for frustum culling. + * + * TODO: In the future, we can add a fast path in the case that the source + * tiling is already in EPSG:3857. */ getBoundingVolume( zRange: ZRange, project: ((xyz: number[]) => number[]) | null, ) { - const tileMatrix = this.tileMatrix; - const { tileWidth, tileHeight, geotransform } = tileMatrix; - - const projectedBounds = computeProjectedTileBounds({ - x: this.x, - y: this.y, - transform: geotransform, - tileWidth, - tileHeight, - }); - - // I have: - // - a tile index x, y, z - // - tileMatrix is at the current z - // - // I need first the projected bounds - - // Sample reference points across the tile surface - const refPoints = REF_POINTS_9; - - console.log("refPoints", refPoints); - - /** Reference points positions in image CRS */ - const refPointPositionsImage: [number, number][] = []; - - for (const [pX, pY] of refPoints) { - // pX, pY are in [0, 1] range - // Interpolate pixel coordinates within the tile - const col = pixelMinCol + pX * (pixelMaxCol - pixelMinCol); - const row = pixelMinRow + pY * (pixelMaxRow - pixelMinRow); - - // Convert pixel coordinates to geographic coordinates using geotransform - const geoX = a * col + b * row + c; - const geoY = d * col + e * row + f; - - refPointPositionsImage.push([geoX, geoY]); - } - - console.log("refPointPositionsImage (image CRS):", refPointPositionsImage); - console.log("Geotransform [a,b,c,d,e,f]:", [a, b, c, d, e, f]); - + // Case 1: Globe view - need to construct an oriented bounding box from + // reprojected sample points, but also using the `project` param if (project) { - assert( - false, - "TODO: implement bounding volume implementation in Globe view", - ); + assert(false, "TODO: implement getBoundingVolume in Globe view"); // Reproject positions to wgs84 instead, then pass them into `project` // return makeOrientedBoundingBoxFromPoints(refPointPositions); } + // (Future) Case 2: Web Mercator input image, can directly compute AABB in + // common space + + // (Future) Case 3: Source projection is already mercator, like UTM. We + // don't need to sample from reference points, we can only use the 4 + // corners. + + // Case 4: Generic case - sample reference points and reproject to + // Web Mercator, then convert to deck.gl common space + return this._getGenericBoundingVolume(zRange); + /** Reference points positions in EPSG 3857 */ const refPointPositionsProjected: [number, number][] = []; @@ -462,6 +435,30 @@ export class RasterTileNode { return obb; } + + /** + * Generic case - sample reference points and reproject to Web Mercator, then + * convert to deck.gl common space + * + */ + _getGenericBoundingVolume(zRange: ZRange): OrientedBoundingBox { + const tileMatrix = this.tileMatrix; + const { tileWidth, tileHeight, geotransform } = tileMatrix; + + const projectedBounds = computeProjectedTileBounds({ + x: this.x, + y: this.y, + transform: geotransform, + tileWidth, + tileHeight, + }); + + const refPointsEPSG3857 = sampleReferencePointsInEPSG3857( + REF_POINTS_9, + projectedBounds, + this.metadata.projectTo3857, + ); + } } /** @@ -512,7 +509,53 @@ function computeProjectedTileBounds({ const maxX = a * pixelMaxCol + b * pixelMaxRow + c; const maxY = d * pixelMaxCol + e * pixelMaxRow + f; - return [minX, minY, maxX, maxY]; + // Note: often `e` in the geotransform is negative (for a north up image when + // the origin is in the **top** left, then increasing the pixel row means + // going down in geospatial space), so maxY < minY + // + // We want to always return an axis-aligned bbox in the form of + // [minX, minY, maxX, maxY], so we need to swap if necessary. + // + // For now, we just use Math.min/Math.max to ensure correct ordering, but we + // could remove the min/max calls if we assume that `a` and `e` are always + // positive/negative respectively. + return [ + Math.min(minX, maxX), + Math.min(minY, maxY), + Math.max(minX, maxX), + Math.max(minY, maxY), + ]; +} + +/** + * Sample the selected reference points in EPSG:3857 + * + * Note that EPSG:3857 is **not** the same as deck.gl's common space! deck.gl's + * common space is the size of `TILE_SIZE` (512) units, while EPSG:3857 uses + * meters. + * + * @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] + */ +function sampleReferencePointsInEPSG3857( + refPoints: [number, number][], + tileBounds: [number, number, number, number], + projectTo3857: (xy: [number, number]) => [number, number], +): [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); + + // Reproject to Web Mercator (EPSG 3857) + const projected = projectTo3857([geoX, geoY]); + refPointPositions.push(projected); + } + + return refPointPositions; +} } /** From ba9cbca1f8e9155d4eb4fc8f06e07c5ff6d27b97 Mon Sep 17 00:00:00 2001 From: Kyle Barron Date: Tue, 16 Dec 2025 14:37:34 -0500 Subject: [PATCH 13/14] implement bounding volume --- .../raster-tileset/raster-tile-traversal.ts | 149 +++++++++++------- 1 file changed, 90 insertions(+), 59 deletions(-) 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 9f9a347..6e8e515 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 @@ -82,6 +82,18 @@ const REF_POINTS_9 = REF_POINTS_5.concat([ [0.5, 1], // bottom edge ]); +/** semi-major axis of the WGS84 ellipsoid + * + * EPSG:3857 also uses the WGS84 datum, so this is used for conversions from + * 3857 to deck.gl common space (scaled to [0-512]) + */ +const WGS84_ELLIPSOID_A = 6378137; + +/** + * Full circumference of the EPSG:3857 Web Mercator world, in meters + */ +const EPSG_3857_CIRCUMFERENCE = 2 * Math.PI * WGS84_ELLIPSOID_A; + /** * Raster Tile Node - represents a single tile in the TileMatrixSet structure * @@ -367,73 +379,56 @@ export class RasterTileNode { // Web Mercator, then convert to deck.gl common space return this._getGenericBoundingVolume(zRange); - /** Reference points positions in EPSG 3857 */ - const refPointPositionsProjected: [number, number][] = []; + // /** Reference points positions in EPSG 3857 */ + // const refPointPositionsProjected: [number, number][] = []; - for (const [pX, pY] of refPointPositionsImage) { - // Reproject to Web Mercator (EPSG 3857) - const projected = this.metadata.projectTo3857([pX, pY]); - refPointPositionsProjected.push(projected); + // // Convert from Web Mercator meters to deck.gl's common space (world units) + // // Web Mercator range: [-20037508.34, 20037508.34] meters + // // deck.gl world space: [0, 512] + // const WEB_MERCATOR_MAX = 20037508.342789244; // Half Earth circumference - // Also log WGS84 for comparison - const wgs84 = this.metadata.projectToWgs84([pX, pY]); - console.log( - `Image [${pX.toFixed(2)}, ${pY.toFixed(2)}] -> WGS84 [${wgs84[0].toFixed(6)}, ${wgs84[1].toFixed(6)}] -> WebMerc [${projected[0].toFixed(2)}, ${projected[1].toFixed(2)}]`, - ); - } + // /** Reference points positions in deck.gl world space */ + // const refPointPositionsWorld: [number, number, number][] = []; - console.log( - "refPointPositionsProjected (EPSG:3857):", - refPointPositionsProjected, - ); - - // Convert from Web Mercator meters to deck.gl's common space (world units) - // Web Mercator range: [-20037508.34, 20037508.34] meters - // deck.gl world space: [0, 512] - const WEB_MERCATOR_MAX = 20037508.342789244; // Half Earth circumference - - /** Reference points positions in deck.gl world space */ - const refPointPositionsWorld: [number, number, number][] = []; + // for (const [mercX, mercY] of refPointPositionsProjected) { + // // X: offset from [-20M, 20M] to [0, 40M], then normalize to [0, 512] + // const worldX = + // ((mercX + WEB_MERCATOR_MAX) / (2 * WEB_MERCATOR_MAX)) * TILE_SIZE; - for (const [mercX, mercY] of refPointPositionsProjected) { - // X: offset from [-20M, 20M] to [0, 40M], then normalize to [0, 512] - const worldX = - ((mercX + WEB_MERCATOR_MAX) / (2 * WEB_MERCATOR_MAX)) * TILE_SIZE; + // // Y: same transformation WITHOUT flip + // // Testing hypothesis: Y-flip might be incorrect since geotransform already handles orientation + // const worldY = + // ((mercY + WEB_MERCATOR_MAX) / (2 * WEB_MERCATOR_MAX)) * TILE_SIZE; - // Y: same transformation WITHOUT flip - // Testing hypothesis: Y-flip might be incorrect since geotransform already handles orientation - const worldY = - ((mercY + WEB_MERCATOR_MAX) / (2 * WEB_MERCATOR_MAX)) * TILE_SIZE; + // console.log( + // `WebMerc [${mercX.toFixed(2)}, ${mercY.toFixed(2)}] -> World [${worldX.toFixed(4)}, ${worldY.toFixed(4)}]`, + // ); - console.log( - `WebMerc [${mercX.toFixed(2)}, ${mercY.toFixed(2)}] -> World [${worldX.toFixed(4)}, ${worldY.toFixed(4)}]`, - ); - - // Add z-range minimum - refPointPositionsWorld.push([worldX, worldY, zRange[0]]); - } - - // Add top z-range if elevation varies - if (zRange[0] !== zRange[1]) { - for (const [mercX, mercY] of refPointPositionsProjected) { - const worldX = - ((mercX + WEB_MERCATOR_MAX) / (2 * WEB_MERCATOR_MAX)) * TILE_SIZE; - const worldY = - TILE_SIZE - - ((mercY + WEB_MERCATOR_MAX) / (2 * WEB_MERCATOR_MAX)) * TILE_SIZE; + // // Add z-range minimum + // refPointPositionsWorld.push([worldX, worldY, zRange[0]]); + // } - refPointPositionsWorld.push([worldX, worldY, zRange[1]]); - } - } + // // Add top z-range if elevation varies + // if (zRange[0] !== zRange[1]) { + // for (const [mercX, mercY] of refPointPositionsProjected) { + // const worldX = + // ((mercX + WEB_MERCATOR_MAX) / (2 * WEB_MERCATOR_MAX)) * TILE_SIZE; + // const worldY = + // TILE_SIZE - + // ((mercY + WEB_MERCATOR_MAX) / (2 * WEB_MERCATOR_MAX)) * TILE_SIZE; + + // refPointPositionsWorld.push([worldX, worldY, zRange[1]]); + // } + // } - console.log("refPointPositionsWorld", refPointPositionsWorld); - console.log("zRange used:", zRange); + // console.log("refPointPositionsWorld", refPointPositionsWorld); + // console.log("zRange used:", zRange); - const obb = makeOrientedBoundingBoxFromPoints(refPointPositionsWorld); - console.log("Created OBB center:", obb.center); - console.log("Created OBB halfAxes:", obb.halfAxes); + // const obb = makeOrientedBoundingBoxFromPoints(refPointPositionsWorld); + // console.log("Created OBB center:", obb.center); + // console.log("Created OBB halfAxes:", obb.halfAxes); - return obb; + // return obb; } /** @@ -444,8 +439,9 @@ export class RasterTileNode { _getGenericBoundingVolume(zRange: ZRange): OrientedBoundingBox { const tileMatrix = this.tileMatrix; const { tileWidth, tileHeight, geotransform } = tileMatrix; + const [minZ, maxZ] = zRange; - const projectedBounds = computeProjectedTileBounds({ + const tileCrsBounds = computeProjectedTileBounds({ x: this.x, y: this.y, transform: geotransform, @@ -455,9 +451,25 @@ export class RasterTileNode { const refPointsEPSG3857 = sampleReferencePointsInEPSG3857( REF_POINTS_9, - projectedBounds, + tileCrsBounds, this.metadata.projectTo3857, ); + + const commonSpacePositions = refPointsEPSG3857.map((xy) => + rescaleEPSG3857ToCommonSpace(xy), + ); + + const refPointPositions: [number, number, number][] = []; + for (const p of commonSpacePositions) { + refPointPositions.push([p[0], p[1], minZ]); + + if (minZ !== maxZ) { + // Also sample at maximum elevation to capture the full 3D volume + refPointPositions.push([p[0], p[1], maxZ]); + } + } + + return makeOrientedBoundingBoxFromPoints(refPointPositions); } } @@ -556,6 +568,25 @@ function sampleReferencePointsInEPSG3857( return refPointPositions; } + +/** + * Rescale positions from EPSG:3857 into deck.gl's common space + * + * Similar to the upstream code here: + * https://github.com/visgl/deck.gl/blob/b0134f025148b52b91320d16768ab5d14a745328/modules/geo-layers/src/tileset-2d/tile-2d-traversal.ts#L172-L177 + * + * @param {number[]} xy [xy description] + * + * @return {number} [return description] + */ +function rescaleEPSG3857ToCommonSpace([x, y]: [number, number]): [ + number, + number, +] { + return [ + (x / EPSG_3857_CIRCUMFERENCE + 0.5) * TILE_SIZE, + (0.5 - y / EPSG_3857_CIRCUMFERENCE) * TILE_SIZE, + ]; } /** From 1e169217fcf0a10489a2fe6784478e02e28f51cd Mon Sep 17 00:00:00 2001 From: Kyle Barron Date: Tue, 16 Dec 2025 14:37:54 -0500 Subject: [PATCH 14/14] clamp y to web mercator bounds when converting epsg:3857 to common space --- .../src/raster-tileset/raster-tile-traversal.ts | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) 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 6e8e515..51479a3 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 @@ -93,6 +93,7 @@ const WGS84_ELLIPSOID_A = 6378137; * Full circumference of the EPSG:3857 Web Mercator world, in meters */ const EPSG_3857_CIRCUMFERENCE = 2 * Math.PI * WGS84_ELLIPSOID_A; +const EPSG_3857_HALF_CIRCUMFERENCE = EPSG_3857_CIRCUMFERENCE / 2; /** * Raster Tile Node - represents a single tile in the TileMatrixSet structure @@ -583,9 +584,15 @@ function rescaleEPSG3857ToCommonSpace([x, y]: [number, number]): [ number, number, ] { + // Clamp Y to Web Mercator bounds + const clampedY = Math.max( + -EPSG_3857_HALF_CIRCUMFERENCE, + Math.min(EPSG_3857_HALF_CIRCUMFERENCE, y), + ); + return [ (x / EPSG_3857_CIRCUMFERENCE + 0.5) * TILE_SIZE, - (0.5 - y / EPSG_3857_CIRCUMFERENCE) * TILE_SIZE, + (0.5 - clampedY / EPSG_3857_CIRCUMFERENCE) * TILE_SIZE, ]; }