diff --git a/examples/cog-basic/package.json b/examples/cog-basic/package.json index 63768fc..c1f2cba 100644 --- a/examples/cog-basic/package.json +++ b/examples/cog-basic/package.json @@ -14,12 +14,11 @@ "@deck.gl/layers": "^9.2.7", "@deck.gl/mapbox": "^9.2.7", "@deck.gl/mesh-layers": "^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", - "geotiff": "^2.1.3", - "geotiff-geokeys-to-proj4": "^2024.4.13", "maplibre-gl": "^5.17.0", "proj4": "^2.20.2", "react": "^19.2.4", diff --git a/examples/cog-basic/src/App.tsx b/examples/cog-basic/src/App.tsx index c0cabb1..1727adb 100644 --- a/examples/cog-basic/src/App.tsx +++ b/examples/cog-basic/src/App.tsx @@ -1,7 +1,6 @@ import type { DeckProps } from "@deck.gl/core"; import { MapboxOverlay } from "@deck.gl/mapbox"; -import { COGLayer, proj } from "@developmentseed/deck.gl-geotiff"; -import { toProj4 } from "geotiff-geokeys-to-proj4"; +import { COGLayer } from "@developmentseed/deck.gl-geotiff"; import "maplibre-gl/dist/maplibre-gl.css"; import { useRef, useState } from "react"; import type { MapRef } from "react-map-gl/maplibre"; @@ -13,26 +12,14 @@ function DeckGLOverlay(props: DeckProps) { return null; } -async function geoKeysParser( - geoKeys: Record, -): Promise { - const projDefinition = toProj4(geoKeys as any); - - return { - def: projDefinition.proj4, - parsed: proj.parseCrs(projDefinition.proj4), - coordinatesUnits: projDefinition.coordinatesUnits as proj.SupportedCrsUnit, - }; -} - -// const COG_URL = -// "https://nz-imagery.s3-ap-southeast-2.amazonaws.com/new-zealand/new-zealand_2024-2025_10m/rgb/2193/CC11.tiff"; +const COG_URL = + "https://nz-imagery.s3-ap-southeast-2.amazonaws.com/new-zealand/new-zealand_2024-2025_10m/rgb/2193/CC11.tiff"; // const COG_URL = // "https://ds-wheels.s3.us-east-1.amazonaws.com/m_4007307_sw_18_060_20220803.tif"; -const COG_URL = - "https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/18/T/WL/2026/1/S2B_18TWL_20260101_0_L2A/TCI.tif"; +// const COG_URL = +// "https://sentinel-cogs.s3.us-west-2.amazonaws.com/sentinel-s2-l2a-cogs/18/T/WL/2026/1/S2B_18TWL_20260101_0_L2A/TCI.tif"; // const COG_URL = // "https://ds-wheels.s3.us-east-1.amazonaws.com/Annual_NLCD_LndCov_2023_CU_C1V0.tif"; @@ -62,7 +49,6 @@ export default function App() { geotiff: COG_URL, debug, debugOpacity, - geoKeysParser, onGeoTIFFLoad: (tiff, options) => { (window as any).tiff = tiff; const { west, south, east, north } = options.geographicBounds; diff --git a/packages/affine/src/affine.ts b/packages/affine/src/affine.ts index b063ca2..53b96a1 100644 --- a/packages/affine/src/affine.ts +++ b/packages/affine/src/affine.ts @@ -81,7 +81,11 @@ export function invert([sa, sb, sc, sd, se, sf]: Affine): Affine { const rd = -sd * idet; const re = sa * idet; - return [ra, rb, -sc * ra - sf * rb, rd, re, -sc * rd - sf * re]; + // biome-ignore format: array + return [ + ra, rb, -sc * ra - sf * rb, + rd, re, -sc * rd - sf * re + ]; } /** Get the 'a' component of an Affine transform. */ diff --git a/packages/deck.gl-geotiff/package.json b/packages/deck.gl-geotiff/package.json index cba26f6..3378363 100644 --- a/packages/deck.gl-geotiff/package.json +++ b/packages/deck.gl-geotiff/package.json @@ -46,11 +46,15 @@ "vitest": "^4.0.18" }, "dependencies": { + "@cogeotiff/core": "^9.1.2", + "@developmentseed/affine": "workspace:^", "@developmentseed/deck.gl-raster": "workspace:^", + "@developmentseed/geotiff": "workspace:^", + "@developmentseed/morecantile": "workspace:^", "@developmentseed/raster-reproject": "workspace:^", "flatbush": "^4.5.0", - "geotiff": "2.1.3", - "proj4": "^2.20.2" + "proj4": "^2.20.2", + "wkt-parser": "^1.5.2" }, "peerDependencies": { "@deck.gl/core": "^9.2.7", diff --git a/packages/deck.gl-geotiff/src/cog-layer.ts b/packages/deck.gl-geotiff/src/cog-layer.ts index 0e481c8..4fb325e 100644 --- a/packages/deck.gl-geotiff/src/cog-layer.ts +++ b/packages/deck.gl-geotiff/src/cog-layer.ts @@ -13,27 +13,26 @@ import type { } from "@deck.gl/geo-layers"; import { TileLayer } from "@deck.gl/geo-layers"; import { PathLayer } from "@deck.gl/layers"; -import type { - RasterModule, - TileMatrix, - TileMatrixSet, +import type { RasterModule } from "@developmentseed/deck.gl-raster"; +import { + RasterLayer, + TileMatrixSetTileset, } from "@developmentseed/deck.gl-raster"; -import { RasterLayer, RasterTileset2D } from "@developmentseed/deck.gl-raster"; +import type { GeoTIFF, Overview } from "@developmentseed/geotiff"; +import { generateTileMatrixSet } from "@developmentseed/geotiff"; +import type { TileMatrixSet } from "@developmentseed/morecantile"; +import { tileTransform } from "@developmentseed/morecantile"; import type { ReprojectionFns } from "@developmentseed/raster-reproject"; import type { Device } from "@luma.gl/core"; -import type { BaseClient, GeoTIFF, GeoTIFFImage, Pool } from "geotiff"; import proj4 from "proj4"; -import { parseCOGTileMatrixSet } from "./cog-tile-matrix-set.js"; -import { - defaultPool, - fetchGeoTIFF, - getGeographicBounds, -} from "./geotiff/geotiff.js"; +import type { ProjectionDefinition } from "wkt-parser"; +import wktParser from "wkt-parser"; +import { fetchGeoTIFF, getGeographicBounds } from "./geotiff/geotiff.js"; import type { TextureDataT } from "./geotiff/render-pipeline.js"; import { inferRenderPipeline } from "./geotiff/render-pipeline.js"; -import { fromGeoTransform } from "./geotiff-reprojection.js"; -import type { GeoKeysParser, ProjectionInfo } from "./proj.js"; -import { epsgIoGeoKeyParser } from "./proj.js"; +import { fromAffine } from "./geotiff-reprojection.js"; +import type { EpsgResolver } from "./proj.js"; +import { epsgResolver } from "./proj.js"; /** * Minimum interface that **must** be returned from getTileData. @@ -52,14 +51,18 @@ export type GetTileDataOptions = { /** The luma.gl Device */ device: Device; - /** the subset to read data from in pixels. */ - window?: [number, number, number, number]; + /** The x coordinate of the tile within the IFD. */ + x: number; + + /** The y coordinate of the tile within the IFD. */ + y: number; /** An AbortSignal that may be signalled if the request is to be aborted */ signal?: AbortSignal; /** The decoder pool to use. */ - pool: Pool; + // TODO: restore pool with new GeoTIFF backend + // pool: Pool; }; type GetTileDataResult = { @@ -79,17 +82,18 @@ export interface COGLayerProps * - An instance of GeoTIFF.js's GeoTIFF class * - An instance of GeoTIFF.js's BaseClient for custom fetching */ - geotiff: GeoTIFF | string | ArrayBuffer | Blob | BaseClient; + geotiff: GeoTIFF | string | URL | ArrayBuffer; /** - * A function callback for parsing GeoTIFF geo keys to a Proj4 compatible - * definition. + * A function callback for parsing numeric EPSG codes to projection + * information (as returned by `wkt-parser`). * - * By default, uses epsg.io to resolve EPSG codes found in the GeoTIFF. - * Alternatively, you may want to use `geotiff-geokeys-to-proj4`, which is - * more extensive but adds 1.5MB to your bundle size. + * The default implementation: + * - makes a request to epsg.io to resolve EPSG codes found in the GeoTIFF. + * - caches any previous requests + * - parses PROJJSON response with `wkt-parser` */ - geoKeysParser?: GeoKeysParser; + epsgResolver?: EpsgResolver; /** * GeoTIFF.js Pool for decoding image chunks. @@ -97,7 +101,8 @@ export interface COGLayerProps * If none is provided, a default Pool will be created and shared between all * COGLayer and GeoTIFFLayer instances. */ - pool?: Pool; + // TODO: Restore a sort of worker pool with our geotiff implementation + // pool?: Pool; /** * Maximum reprojection error in pixels for mesh refinement. @@ -116,7 +121,7 @@ export interface COGLayerProps * luma.gl, along with optional custom shaders for the RasterLayer. */ getTileData?: ( - image: GeoTIFFImage, + image: GeoTIFF | Overview, options: GetTileDataOptions, ) => Promise; @@ -145,14 +150,11 @@ export interface COGLayerProps /** * Called when the GeoTIFF metadata has been loaded and parsed. - * - * @param {GeoTIFF} geotiff - * @param {ProjectionInfo} projection */ onGeoTIFFLoad?: ( geotiff: GeoTIFF, options: { - projection: ProjectionInfo; + projection: ProjectionDefinition; /** * Bounds of the image in geographic coordinates (WGS84) [minLon, minLat, * maxLon, maxLat] @@ -176,7 +178,7 @@ export interface COGLayerProps } const defaultProps: Partial = { - geoKeysParser: epsgIoGeoKeyParser, + epsgResolver, debug: false, debugOpacity: 0.5, }; @@ -191,10 +193,11 @@ export class COGLayer< static override defaultProps = defaultProps; declare state: { - forwardReproject?: ReprojectionFns["forwardReproject"]; - inverseReproject?: ReprojectionFns["inverseReproject"]; - metadata?: TileMatrixSet; - images?: GeoTIFFImage[]; + geotiff: GeoTIFF; + forwardTo4326?: ReprojectionFns["forwardReproject"]; + inverseFrom4326?: ReprojectionFns["inverseReproject"]; + forwardTo3857?: ReprojectionFns["forwardReproject"]; + tms?: TileMatrixSet; defaultGetTileData?: COGLayerProps["getTileData"]; defaultRenderTile?: COGLayerProps["renderTile"]; }; @@ -218,32 +221,30 @@ export class COGLayer< async _parseGeoTIFF(): Promise { const geotiff = await fetchGeoTIFF(this.props.geotiff); - - const geoKeysParser = this.props.geoKeysParser!; - const metadata = await parseCOGTileMatrixSet(geotiff, geoKeysParser); - - const image = await geotiff.getImage(); - const imageCount = await geotiff.getImageCount(); - const images: GeoTIFFImage[] = []; - for (let imageIdx = 0; imageIdx < imageCount; imageIdx++) { - images.push(await geotiff.getImage(imageIdx)); - } - - const sourceProjection = await geoKeysParser(image.getGeoKeys()); - if (!sourceProjection) { - throw new Error( - "Could not determine source projection from GeoTIFF geo keys", - ); - } - - const converter = proj4(sourceProjection.def, "EPSG:4326"); - const forwardReproject = (x: number, y: number) => - converter.forward<[number, number]>([x, y], false); - const inverseReproject = (x: number, y: number) => - converter.inverse<[number, number]>([x, y], false); + const crs = geotiff.crs; + const sourceProjection = + typeof crs === "number" + ? await this.props.epsgResolver!(crs) + : wktParser(crs); + + const tms = generateTileMatrixSet(geotiff, sourceProjection); + + // @ts-expect-error - proj4 typings are incomplete and don't support + // wkt-parser input + const converter4326 = proj4(sourceProjection, "EPSG:4326"); + const forwardTo4326 = (x: number, y: number) => + converter4326.forward<[number, number]>([x, y], false); + const inverseFrom4326 = (x: number, y: number) => + converter4326.inverse<[number, number]>([x, y], false); + + // @ts-expect-error - proj4 typings are incomplete and don't support + // wkt-parser input + const converter3857 = proj4(sourceProjection, "EPSG:3857"); + const forwardTo3857 = (x: number, y: number) => + converter3857.forward<[number, number]>([x, y], false); if (this.props.onGeoTIFFLoad) { - const geographicBounds = getGeographicBounds(image, converter); + const geographicBounds = getGeographicBounds(geotiff, converter4326); this.props.onGeoTIFFLoad(geotiff, { projection: sourceProjection, geographicBounds, @@ -251,13 +252,14 @@ export class COGLayer< } const { getTileData: defaultGetTileData, renderTile: defaultRenderTile } = - inferRenderPipeline(image.fileDirectory, this.context.device); + inferRenderPipeline(geotiff, this.context.device); this.setState({ - metadata, - forwardReproject, - inverseReproject, - images, + geotiff, + tms, + forwardTo4326, + inverseFrom4326, + forwardTo3857, defaultGetTileData, defaultRenderTile, }); @@ -268,30 +270,27 @@ export class COGLayer< */ async _getTileData( tile: TileLoadProps, - images: GeoTIFFImage[], - metadata: TileMatrixSet, + geotiff: GeoTIFF, + tms: TileMatrixSet, ): Promise> { const { signal } = tile; const { x, y, z } = tile.index; // Select overview image - const geotiffImage = images[images.length - 1 - z]!; - const imageHeight = geotiffImage.getHeight(); - const imageWidth = geotiffImage.getWidth(); + // If z=0, use the coarsest overview (which is the last in the array) + // If z=max, use the full-resolution image (which is the first in the array) - const tileMatrix = metadata.tileMatrices[z]!; - const { tileWidth, tileHeight } = tileMatrix; + // TODO: should be able to (micro) optimize this to not create the array + // Something like: + // const image = z === geotiff.overviews.length - 1 ? geotiff : + // geotiff.overviews[geotiff.overviews.length - 1 - z]!; + const images = [geotiff, ...geotiff.overviews]; + const image = images[images.length - 1 - z]!; - const tileGeotransform = computeTileGeotransform(x, y, tileMatrix); - const { forwardTransform, inverseTransform } = - fromGeoTransform(tileGeotransform); + const tileMatrix = tms.tileMatrices[z]!; - const window: [number, number, number, number] = [ - x * tileWidth, - y * tileHeight, - Math.min((x + 1) * tileWidth, imageWidth), - Math.min((y + 1) * tileHeight, imageHeight), - ]; + const tileAffine = tileTransform(tileMatrix, { col: x, row: y }); + const { forwardTransform, inverseTransform } = fromAffine(tileAffine); const getTileData = this.props.getTileData || this.state.defaultGetTileData!; @@ -302,11 +301,13 @@ export class COGLayer< ? AbortSignal.any([signal, this.props.signal]) : signal || this.props.signal; - const data = await getTileData(geotiffImage, { + const data = await getTileData(image, { device: this.context.device, - window, + x, + y, signal: combinedSignal, - pool: this.props.pool || defaultPool(), + // TODO: restore pool + // pool: this.props.pool || defaultPool(), }); return { @@ -326,9 +327,9 @@ export class COGLayer< _offset: number; tile: Tile2DHeader>; }, - metadata: TileMatrixSet, - forwardReproject: ReprojectionFns["forwardReproject"], - inverseReproject: ReprojectionFns["inverseReproject"], + tms: TileMatrixSet, + forwardTo4326: ReprojectionFns["forwardReproject"], + inverseFrom4326: ReprojectionFns["inverseReproject"], ): Layer | LayersList | null { const { maxError, debug, debugOpacity } = this.props; const { tile } = props; @@ -355,8 +356,8 @@ export class COGLayer< reprojectionFns: { forwardTransform, inverseTransform, - forwardReproject, - inverseReproject, + forwardReproject: forwardTo4326, + inverseReproject: inverseFrom4326, }, debug, debugOpacity, @@ -368,19 +369,24 @@ export class COGLayer< // 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) { + const projectedBounds: { + topLeft: [number, number]; + topRight: [number, number]; + bottomLeft: [number, number]; + bottomRight: [number, number]; + } = (tile as any)?.projectedBounds; + + if (!projectedBounds || !tms) { 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); + const topLeftWgs84 = forwardTo4326(topLeft[0], topLeft[1]); + const topRightWgs84 = forwardTo4326(topRight[0], topRight[1]); + const bottomRightWgs84 = forwardTo4326(bottomRight[0], bottomRight[1]); + const bottomLeftWgs84 = forwardTo4326(bottomLeft[0], bottomLeft[1]); // Create a closed path around the tile bounds const path = [ @@ -407,37 +413,44 @@ export class COGLayer< return layers; } + /** Define the underlying deck.gl TileLayer. */ renderTileLayer( - metadata: TileMatrixSet, - forwardReproject: ReprojectionFns["forwardReproject"], - inverseReproject: ReprojectionFns["inverseReproject"], - images: GeoTIFFImage[], + tms: TileMatrixSet, + forwardTo4326: ReprojectionFns["forwardReproject"], + inverseFrom4326: ReprojectionFns["inverseReproject"], + forwardTo3857: ReprojectionFns["forwardReproject"], + geotiff: GeoTIFF, ): TileLayer { // Create a factory class that wraps COGTileset2D with the metadata - class RasterTileset2DFactory extends RasterTileset2D { + class TileMatrixSetTilesetFactory extends TileMatrixSetTileset { constructor(opts: Tileset2DProps) { - super(metadata, opts); + super(opts, tms, { + projectTo4326: forwardTo4326, + projectTo3857: forwardTo3857, + }); } } return new TileLayer>({ id: `cog-tile-layer-${this.id}`, - TilesetClass: RasterTileset2DFactory, - getTileData: async (tile) => this._getTileData(tile, images, metadata), + TilesetClass: TileMatrixSetTilesetFactory, + getTileData: async (tile) => this._getTileData(tile, geotiff, tms), renderSubLayers: (props) => - this._renderSubLayers( - props, - metadata, - forwardReproject, - inverseReproject, - ), + this._renderSubLayers(props, tms, forwardTo4326, inverseFrom4326), }); } renderLayers() { - const { forwardReproject, inverseReproject, metadata, images } = this.state; - - if (!forwardReproject || !inverseReproject || !metadata || !images) { + const { forwardTo4326, inverseFrom4326, forwardTo3857, tms, geotiff } = + this.state; + + if ( + !forwardTo4326 || + !inverseFrom4326 || + !forwardTo3857 || + !tms || + !geotiff + ) { return null; } @@ -445,34 +458,11 @@ export class COGLayer< // nullable in any part of function scope, the tileset factory wrapper gives // a type error return this.renderTileLayer( - metadata, - forwardReproject, - inverseReproject, - images, + tms, + forwardTo4326, + inverseFrom4326, + forwardTo3857, + geotiff, ); } } - -/** - * Compute the affine geotransform for this tile. - * - * We need to offset the geotransform for the matrix level by the tile's pixel - * origin. - */ -function computeTileGeotransform( - x: number, - y: number, - tileMatrix: TileMatrix, -): [number, number, number, number, number, number] { - const { tileWidth, tileHeight } = tileMatrix; - - const xPixelOrigin = x * tileWidth; - const yPixelOrigin = y * tileHeight; - - const [a, b, c, d, e, f] = tileMatrix.geotransform; - - const xCoordOffset = a * xPixelOrigin + b * yPixelOrigin + c; - const yCoordOffset = d * xPixelOrigin + e * yPixelOrigin + f; - - return [a, b, xCoordOffset, d, e, yCoordOffset]; -} diff --git a/packages/deck.gl-geotiff/src/cog-tile-matrix-set.ts b/packages/deck.gl-geotiff/src/cog-tile-matrix-set.ts deleted file mode 100644 index 20433cc..0000000 --- a/packages/deck.gl-geotiff/src/cog-tile-matrix-set.ts +++ /dev/null @@ -1,266 +0,0 @@ -import type { - TileMatrix, - TileMatrixSet, - TileMatrixSetBoundingBox, -} from "@developmentseed/deck.gl-raster"; -import type { GeoTIFF, GeoTIFFImage } from "geotiff"; -import proj4, { type ProjectionDefinition } from "proj4"; -import Ellipsoid from "./ellipsoids.js"; -import { extractGeotransform } from "./geotiff-reprojection"; -import type { GeoKeysParser, ProjectionInfo, SupportedCrsUnit } from "./proj"; - -// 0.28 mm per pixel -// https://docs.ogc.org/is/17-083r4/17-083r4.html#toc15 -const SCREEN_PIXEL_SIZE = 0.00028; - -/** - * - * Ported from Vincent's work here: - * https://github.com/developmentseed/morecantile/pull/187/changes#diff-402eedddfa30af554d03750c8a82a09962b44b044976c321b774b484b98e8f48R665 - * - * @return {TileMatrixSet}[return description] - */ -export async function parseCOGTileMatrixSet( - tiff: GeoTIFF, - geoKeysParser: GeoKeysParser, -): Promise { - const fullResImage = await tiff.getImage(); - - if (!fullResImage.isTiled) { - throw new Error("COG TileMatrixSet requires a tiled GeoTIFF"); - } - - const imageCount = await tiff.getImageCount(); - const bbox = fullResImage.getBoundingBox(); - const fullImageWidth = fullResImage.getWidth(); - const fullImageHeight = fullResImage.getHeight(); - - const crs = await geoKeysParser(fullResImage.getGeoKeys()); - - if (crs === null) { - throw new Error( - "Could not determine coordinate reference system from GeoTIFF geo keys", - ); - } - - const projectToWgs84 = proj4(crs.def, "EPSG:4326").forward; - const projectTo3857 = proj4(crs.def, "EPSG:3857").forward; - - const boundingBox: TileMatrixSet["boundingBox"] = { - lowerLeft: [bbox[0]!, bbox[1]!], - upperRight: [bbox[2]!, bbox[3]!], - }; - - const transform = extractGeotransform(fullResImage); - - if (transform[1] !== 0 || transform[3] !== 0) { - // TileMatrixSet assumes orthogonal axes - throw new Error( - "COG TileMatrixSet with rotation/skewed geotransform is not supported", - ); - } - - const cellSize = Math.abs(transform[0]); - - const tileWidth = fullResImage.getTileWidth(); - const tileHeight = fullResImage.getTileHeight(); - - const tileMatrices: TileMatrix[] = [ - { - // Set as highest resolution / finest level - id: String(imageCount - 1), - scaleDenominator: - (cellSize * metersPerUnit(crs.parsed, crs.coordinatesUnits)) / - SCREEN_PIXEL_SIZE, - cellSize, - pointOfOrigin: [transform[2], transform[5]], - tileWidth: fullResImage.getTileWidth(), - tileHeight: fullResImage.getTileHeight(), - matrixWidth: Math.ceil(fullImageWidth / tileWidth), - matrixHeight: Math.ceil(fullImageHeight / tileHeight), - geotransform: transform, - }, - ]; - - // Starting from 1 to skip full res image - for (let imageIdx = 1; imageIdx < imageCount; imageIdx++) { - const image = await tiff.getImage(imageIdx); - - if (!image.isTiled) { - throw new Error("COG TileMatrixSet requires a tiled GeoTIFF"); - } - - const tileMatrix = createOverviewTileMatrix({ - id: String(imageCount - 1 - imageIdx), - image, - fullWidth: fullImageWidth, - fullHeight: fullImageHeight, - baseTransform: transform, - crs, - }); - tileMatrices.push(tileMatrix); - } - - // Reverse to have coarsest level first - tileMatrices.reverse(); - - return { - crs, - boundingBox, - wgsBounds: computeWgs84BoundingBox(boundingBox, projectToWgs84), - tileMatrices, - projectToWgs84, - projectTo3857, - }; -} - -/** - * Coefficient to convert the coordinate reference system (CRS) - * units into meters (metersPerUnit). - * - * From note g in http://docs.opengeospatial.org/is/17-083r2/17-083r2.html#table_2: - * - * > If the CRS uses meters as units of measure for the horizontal dimensions, - * > then metersPerUnit=1; if it has degrees, then metersPerUnit=2pa/360 - * > (a is the Earth maximum radius of the ellipsoid). - * - * If `crsUnit` is provided, it takes precedence over the unit defined in the - * CRS. This exists because sometimes the parsed CRS from - * `geotiff-geokeys-to-proj4` doesn't have a unit defined. - */ -// https://github.com/developmentseed/morecantile/blob/7c95a11c491303700d6e33e9c1607f2719584dec/morecantile/utils.py#L67-L90 -function metersPerUnit( - parsedCrs: ProjectionDefinition, - crsUnit?: SupportedCrsUnit, -): number { - const unit = (crsUnit || parsedCrs.units)?.toLowerCase(); - switch (unit) { - case "m": - case "metre": - case "meter": - case "meters": - return 1; - case "foot": - return 0.3048; - case "us survey foot": - return 1200 / 3937; - } - - if (unit === "degree") { - // 2 * π * ellipsoid semi-major-axis / 360 - const { a } = Ellipsoid[parsedCrs.ellps as keyof typeof Ellipsoid]; - return (2 * Math.PI * a) / 360; - } - - throw new Error(`Unsupported CRS units: ${unit}`); -} - -/** - * Create tile matrix for COG overview - */ -function createOverviewTileMatrix({ - id, - image, - fullWidth, - baseTransform, - crs, -}: { - id: string; - image: GeoTIFFImage; - fullWidth: number; - fullHeight: number; - baseTransform: [number, number, number, number, number, number]; - crs: ProjectionInfo; -}): TileMatrix { - const width = image.getWidth(); - const height = image.getHeight(); - - // For now, just use scaleX - // https://github.com/developmentseed/morecantile/pull/187/changes#r2621314673 - const scaleX = fullWidth / width; - - // const scaleY = fullHeight / height; - // if (Math.abs(scaleX - scaleY) > 1e-3) { - // throw new Error("Non-uniform overview scaling detected (X/Y differ)"); - // } - - const scale = scaleX; - - const geotransform: [number, number, number, number, number, number] = [ - baseTransform[0] * scale, - baseTransform[1] * scale, - baseTransform[2], // x origin stays the same - baseTransform[3] * scale, - baseTransform[4] * scale, - baseTransform[5], // y origin stays the same - ]; - const cellSize = Math.abs(geotransform[0]); - - const tileWidth = image.getTileWidth(); - const tileHeight = image.getTileHeight(); - - return { - id, - scaleDenominator: - (cellSize * metersPerUnit(crs.parsed, crs.coordinatesUnits)) / - SCREEN_PIXEL_SIZE, - cellSize, - pointOfOrigin: [geotransform[2], geotransform[5]], - tileWidth, - tileHeight, - matrixWidth: Math.ceil(width / tileWidth), - matrixHeight: Math.ceil(height / tileHeight), - geotransform, - }; -} - -function computeWgs84BoundingBox( - boundingBox: TileMatrixSetBoundingBox, - projectToWgs84: (point: [number, number]) => [number, number], -): TileMatrixSetBoundingBox { - const lowerLeftWgs84 = projectToWgs84(boundingBox.lowerLeft); - const lowerRightWgs84 = projectToWgs84([ - boundingBox.upperRight[0], - boundingBox.lowerLeft[1], - ]); - const upperRightWgs84 = projectToWgs84(boundingBox.upperRight); - const upperLeftWgs84 = projectToWgs84([ - boundingBox.lowerLeft[0], - boundingBox.upperRight[1], - ]); - - // Compute min/max lat/lon - const minLon = Math.min( - lowerLeftWgs84[0], - lowerRightWgs84[0], - upperRightWgs84[0], - upperLeftWgs84[0], - ); - const maxLon = Math.max( - lowerLeftWgs84[0], - lowerRightWgs84[0], - upperRightWgs84[0], - upperLeftWgs84[0], - ); - const minLat = Math.min( - lowerLeftWgs84[1], - lowerRightWgs84[1], - upperRightWgs84[1], - upperLeftWgs84[1], - ); - const maxLat = Math.max( - lowerLeftWgs84[1], - lowerRightWgs84[1], - upperRightWgs84[1], - upperLeftWgs84[1], - ); - - return { - lowerLeft: [minLon, minLat], - upperRight: [maxLon, maxLat], - }; -} - -export const __TEST_EXPORTS = { - metersPerUnit, -}; diff --git a/packages/deck.gl-geotiff/src/ellipsoids.ts b/packages/deck.gl-geotiff/src/ellipsoids.ts deleted file mode 100644 index 48f73f4..0000000 --- a/packages/deck.gl-geotiff/src/ellipsoids.ts +++ /dev/null @@ -1,154 +0,0 @@ -/** - * Vendored and edited from proj4.js. - * - * In the implementation of metersPerUnit while generating a COG TileMatrixSet, - * we need to know the size of the semi major axis in the case the CRS is in - * degrees. - * - * https://github.com/proj4js/proj4js/blob/e90e5fa6872a1ffc40edb161cbeb4bd5e3bd9db5/lib/constants/Ellipsoid.js - */ - -const ellipsoids = { - MERIT: { - a: 6378137, - }, - SGS85: { - a: 6378136, - }, - GRS80: { - a: 6378137, - }, - IAU76: { - a: 6378140, - }, - airy: { - a: 6377563.396, - b: 6356256.91, - }, - APL4: { - a: 6378137, - }, - NWL9D: { - a: 6378145, - }, - mod_airy: { - a: 6377340.189, - b: 6356034.446, - }, - andrae: { - a: 6377104.43, - }, - aust_SA: { - a: 6378160, - }, - GRS67: { - a: 6378160, - }, - bessel: { - a: 6377397.155, - }, - bess_nam: { - a: 6377483.865, - }, - clrk66: { - a: 6378206.4, - b: 6356583.8, - }, - clrk80: { - a: 6378249.145, - }, - clrk80ign: { - a: 6378249.2, - b: 6356515, - }, - clrk58: { - a: 6378293.645208759, - }, - CPM: { - a: 6375738.7, - }, - delmbr: { - a: 6376428, - }, - engelis: { - a: 6378136.05, - }, - evrst30: { - a: 6377276.345, - }, - evrst48: { - a: 6377304.063, - }, - evrst56: { - a: 6377301.243, - }, - evrst69: { - a: 6377295.664, - }, - evrstSS: { - a: 6377298.556, - }, - fschr60: { - a: 6378166, - }, - fschr60m: { - a: 6378155, - }, - fschr68: { - a: 6378150, - }, - helmert: { - a: 6378200, - }, - hough: { - a: 6378270, - }, - intl: { - a: 6378388, - }, - kaula: { - a: 6378163, - }, - lerch: { - a: 6378139, - }, - mprts: { - a: 6397300, - }, - new_intl: { - a: 6378157.5, - b: 6356772.2, - }, - plessis: { - a: 6376523, - }, - krass: { - a: 6378245, - }, - SEasia: { - a: 6378155, - b: 6356773.3205, - }, - walbeck: { - a: 6376896, - b: 6355834.8467, - }, - WGS60: { - a: 6378165, - }, - WGS66: { - a: 6378145, - }, - WGS7: { - a: 6378135, - }, - WGS84: { - a: 6378137, - }, - sphere: { - a: 6370997, - b: 6370997, - }, -}; - -export default ellipsoids; diff --git a/packages/deck.gl-geotiff/src/geotiff-layer.ts b/packages/deck.gl-geotiff/src/geotiff-layer.ts index cf7d7f4..ae8f63b 100644 --- a/packages/deck.gl-geotiff/src/geotiff-layer.ts +++ b/packages/deck.gl-geotiff/src/geotiff-layer.ts @@ -1,18 +1,15 @@ import type { CompositeLayerProps, UpdateParameters } from "@deck.gl/core"; import { CompositeLayer } from "@deck.gl/core"; import { RasterLayer } from "@developmentseed/deck.gl-raster"; +import type { GeoTIFF } from "@developmentseed/geotiff"; import type { ReprojectionFns } from "@developmentseed/raster-reproject"; -import type { BaseClient, GeoTIFF, Pool } from "geotiff"; import proj4 from "proj4"; -import { - defaultPool, - fetchGeoTIFF, - getGeographicBounds, - loadRgbImage, -} from "./geotiff/geotiff.js"; +import type { ProjectionDefinition } from "wkt-parser"; +import wktParser from "wkt-parser"; +import { fetchGeoTIFF, getGeographicBounds } from "./geotiff/geotiff.js"; import { extractGeotiffReprojectors } from "./geotiff-reprojection.js"; -import type { GeoKeysParser, ProjectionInfo } from "./proj.js"; -import { epsgIoGeoKeyParser } from "./proj.js"; +import type { EpsgResolver } from "./proj.js"; +import { epsgResolver } from "./proj.js"; export interface GeoTIFFLayerProps extends CompositeLayerProps { /** @@ -24,17 +21,18 @@ export interface GeoTIFFLayerProps extends CompositeLayerProps { * - An instance of GeoTIFF.js's GeoTIFF class * - An instance of GeoTIFF.js's BaseClient for custom fetching */ - geotiff: GeoTIFF | string | ArrayBuffer | Blob | BaseClient; + geotiff: GeoTIFF | string | URL | ArrayBuffer; /** - * A function callback for parsing GeoTIFF geo keys to a Proj4 compatible - * definition. + * A function callback for parsing numeric EPSG codes to projection + * information (as returned by `wkt-parser`). * - * By default, uses epsg.io to resolve EPSG codes found in the GeoTIFF. - * Alternatively, you may want to use `geotiff-geokeys-to-proj4`, which is - * more extensive but adds 1.5MB to your bundle size. + * The default implementation: + * - makes a request to epsg.io to resolve EPSG codes found in the GeoTIFF. + * - caches any previous requests + * - parses PROJJSON response with `wkt-parser` */ - geoKeysParser?: GeoKeysParser; + epsgResolver?: EpsgResolver; /** * GeoTIFF.js Pool for decoding image chunks. @@ -42,7 +40,8 @@ export interface GeoTIFFLayerProps extends CompositeLayerProps { * If none is provided, a default Pool will be created and shared between all * COGLayer and GeoTIFFLayer instances. */ - pool?: Pool; + // TODO: Restore a sort of worker pool with our geotiff implementation + // pool?: Pool; /** * Maximum reprojection error in pixels for mesh refinement. @@ -65,14 +64,11 @@ export interface GeoTIFFLayerProps extends CompositeLayerProps { /** * Called when the GeoTIFF metadata has been loaded and parsed. - * - * @param {GeoTIFF} geotiff - * @param {ProjectionInfo} projection */ onGeoTIFFLoad?: ( geotiff: GeoTIFF, options: { - projection: ProjectionInfo; + projection: ProjectionDefinition; /** * Bounds of the image in geographic coordinates (WGS84) [minLon, minLat, * maxLon, maxLat] @@ -88,7 +84,7 @@ export interface GeoTIFFLayerProps extends CompositeLayerProps { } const defaultProps = { - geoKeysParser: epsgIoGeoKeyParser, + epsgResolver, }; /** @@ -130,40 +126,42 @@ export class GeoTIFFLayer extends CompositeLayer { async _parseGeoTIFF(): Promise { const geotiff = await fetchGeoTIFF(this.props.geotiff); - const image = await geotiff.getImage(); - - const geoKeysParser = this.props.geoKeysParser!; - const sourceProjection = await geoKeysParser(image.getGeoKeys()); - if (!sourceProjection) { - throw new Error( - "Could not determine source projection from GeoTIFF geo keys", - ); - } + const crs = geotiff.crs; + const sourceProjection = + typeof crs === "number" + ? await this.props.epsgResolver!(crs) + : wktParser(crs); - const converter = proj4(sourceProjection.def, "EPSG:4326"); + // @ts-expect-error proj4 has incomplete types that don't support wkt-parser + // output + const converter = proj4(sourceProjection, "EPSG:4326"); if (this.props.onGeoTIFFLoad) { - const geographicBounds = getGeographicBounds(image, converter); + const geographicBounds = getGeographicBounds(geotiff, converter); this.props.onGeoTIFFLoad(geotiff, { projection: sourceProjection, geographicBounds, }); } + // @ts-expect-error unused variable + // biome-ignore lint/correctness/noUnusedVariables: not implemented const reprojectionFns = await extractGeotiffReprojectors( geotiff, - sourceProjection.def, + sourceProjection, ); - const { texture, height, width } = await loadRgbImage(image, { - pool: this.props.pool || defaultPool(), - }); - - this.setState({ - reprojectionFns, - imageData: texture, - height, - width, - }); + + // Our GeoTIFF implementation doesn't currently support reading the full + // image; it only supports reading tiles. + throw new Error("Loading GeoTIFF image data not yet implemented"); + // const { texture, height, width } = await loadRgbImage(image); + + // this.setState({ + // reprojectionFns, + // imageData: texture, + // height, + // width, + // }); } renderLayers() { diff --git a/packages/deck.gl-geotiff/src/geotiff-reprojection.ts b/packages/deck.gl-geotiff/src/geotiff-reprojection.ts index 379d0b9..457c693 100644 --- a/packages/deck.gl-geotiff/src/geotiff-reprojection.ts +++ b/packages/deck.gl-geotiff/src/geotiff-reprojection.ts @@ -1,106 +1,24 @@ -/* eslint-env browser */ - +import * as affine from "@developmentseed/affine"; +import type { GeoTIFF } from "@developmentseed/geotiff"; import type { ReprojectionFns } from "@developmentseed/raster-reproject"; -import { - applyAffine, - invertGeoTransform, -} from "@developmentseed/raster-reproject/affine"; -import type { GeoTIFF, GeoTIFFImage } from "geotiff"; import proj4 from "proj4"; import type { PROJJSONDefinition } from "proj4/dist/lib/core"; import type Projection from "proj4/dist/lib/Proj"; - -export const OGC_84: PROJJSONDefinition = { - $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, - }, - // @ts-expect-error - proj4 types are incomplete - id: { authority: "OGC", code: "CRS84" }, -}; +import type { ProjectionDefinition } from "wkt-parser"; // Derived from existing work here: // https://github.com/developmentseed/lonboard/blob/35a1f3d691604ad9e083bf10a4bfde4158171486/src/cog-tileset/claude-tileset-2d-improved.ts#L141 // // TODO: return a RasterReprojector instance, given the IFD and tile of interest? export async function extractGeotiffReprojectors( - tiff: GeoTIFF, - sourceProjection: string | PROJJSONDefinition, - outputCrs: string | PROJJSONDefinition | Projection = OGC_84, + geotiff: GeoTIFF, + sourceProjection: string | PROJJSONDefinition | ProjectionDefinition, + outputCrs: string | PROJJSONDefinition | Projection = "EPSG:4326", ): Promise { - const image = await tiff.getImage(); - - // 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); - + // @ts-expect-error - proj4 type definitions are incomplete and don't include + // support for wkt-parser output const converter = proj4(sourceProjection, outputCrs); - const { forwardTransform, inverseTransform } = - fromGeoTransform(baseGeotransform); + const { forwardTransform, inverseTransform } = fromAffine(geotiff.transform); return { forwardTransform, @@ -112,66 +30,17 @@ export async function extractGeotiffReprojectors( }; } -export function fromGeoTransform( +export function fromAffine( geotransform: [number, number, number, number, number, number], ): { forwardTransform: (x: number, y: number) => [number, number]; inverseTransform: (x: number, y: number) => [number, number]; } { - const inverseGeotransform = invertGeoTransform(geotransform); + const inverseGeotransform = affine.invert(geotransform); return { - forwardTransform: (x: number, y: number) => applyAffine(x, y, geotransform), + forwardTransform: (x: number, y: number) => + affine.apply(geotransform, x, y), inverseTransform: (x: number, y: number) => - applyAffine(x, y, inverseGeotransform), + affine.apply(inverseGeotransform, x, y), }; } - -/** - * Extract affine geotransform from a GeoTIFF image. - * - * The first `image` must be passed in, as only the top-level IFD contains geo - * keys. - * - * Returns a 6-element array in Python `affine` package ordering: - * [a, b, c, d, e, f] where: - * - x_geo = a * col + b * row + c - * - y_geo = d * col + e * row + f - * - * This is NOT GDAL ordering, which is [c, a, b, f, d, e]. - */ -export function extractGeotransform( - image: GeoTIFFImage, -): [number, number, number, number, number, number] { - const origin = image.getOrigin(); - const resolution = image.getResolution(); - - // origin: [x, y, z] - // resolution: [x_res, y_res, z_res] - - // Check for rotation/skew in the file directory - const fileDirectory = image.getFileDirectory(); - const modelTransformation = fileDirectory.ModelTransformation; - - let b = 0; // row rotation - let d = 0; // column rotation - - if (modelTransformation && modelTransformation.length >= 16) { - // ModelTransformation is a 4x4 matrix in row-major order - // [0 1 2 3 ] [a b 0 c] - // [4 5 6 7 ] = [d e 0 f] - // [8 9 10 11] [0 0 1 0] - // [12 13 14 15] [0 0 0 1] - b = modelTransformation[1]; - d = modelTransformation[4]; - } - - // Return in affine package ordering: [a, b, c, d, e, f] - return [ - resolution[0]!, // a: pixel width - b, // b: row rotation - origin[0]!, // c: x origin - d, // d: column rotation - resolution[1]!, // e: pixel height (often negative) - origin[1]!, // f: y origin - ]; -} diff --git a/packages/deck.gl-geotiff/src/geotiff/geotiff.ts b/packages/deck.gl-geotiff/src/geotiff/geotiff.ts index 27a1698..58bd347 100644 --- a/packages/deck.gl-geotiff/src/geotiff/geotiff.ts +++ b/packages/deck.gl-geotiff/src/geotiff/geotiff.ts @@ -1,107 +1,49 @@ // Utilities for interacting with geotiff.js. -import type { GeoTIFF, GeoTIFFImage, TypedArrayWithDimensions } from "geotiff"; -import { - BaseClient, - fromArrayBuffer, - fromBlob, - fromCustomClient, - fromUrl, - Pool, -} from "geotiff"; +import type { RasterArray } from "@developmentseed/geotiff"; +import { GeoTIFF } from "@developmentseed/geotiff"; import type { Converter } from "proj4"; -/** - * Options that may be passed when reading image data from geotiff.js - */ -type ReadRasterOptions = { - /** the subset to read data from in pixels. */ - window?: [number, number, number, number]; - - /** The optional decoder pool to use. */ - pool?: Pool; - - /** An AbortSignal that may be signalled if the request is to be aborted */ - signal?: AbortSignal; -}; - -/** - * A default geotiff.js decoder pool instance. - * - * It will be created on first call of `defaultPool`. - */ -let DEFAULT_POOL: Pool | null = null; - -/** - * Retrieve the default geotiff.js decoder Pool. - * - * If a Pool has not yet been created, it will be created on first call. - * - * The Pool will be shared between all COGLayer and GeoTIFFLayer instances. - */ -export function defaultPool(): Pool { - if (DEFAULT_POOL === null) { - DEFAULT_POOL = new Pool(); - } - - return DEFAULT_POOL; -} - -/** - * Load an RGBA image from a GeoTIFFImage. - */ -export async function loadRgbImage( - image: GeoTIFFImage, - options?: ReadRasterOptions, -): Promise<{ texture: ImageData; height: number; width: number }> { - const mergedOptions = { - ...options, - interleave: true, - enableAlpha: true, - }; - // Since we set interleave: true, the result is a single array with all - // samples, so we cast to TypedArrayWithDimensions - // https://github.com/geotiffjs/geotiff.js/issues/486 - const rgbImage = (await image.readRGB( - mergedOptions, - )) as TypedArrayWithDimensions; - const imageData = addAlphaChannel(rgbImage); - - return { - texture: imageData, - height: rgbImage.height, - width: rgbImage.width, - }; -} - /** * Add an alpha channel to an RGB image array. * * Only supports input arrays with 3 (RGB) or 4 (RGBA) channels. If the input is * already RGBA, it is returned unchanged. */ -export function addAlphaChannel(rgbImage: TypedArrayWithDimensions): ImageData { +export function addAlphaChannel(rgbImage: RasterArray): RasterArray { const { height, width } = rgbImage; - if (rgbImage.length === height * width * 4) { + if (rgbImage.layout === "band-separate") { + // This should be pretty easy to do by just returning an additional array of + // 255s + // But not sure if we'll want to do that, because it's fine to upload 3 + // separate textures. + throw new Error("Band-separate images not yet implemented."); + } + + if (rgbImage.data.length === height * width * 4) { // Already has alpha channel - return new ImageData(new Uint8ClampedArray(rgbImage), width, height); - } else if (rgbImage.length === height * width * 3) { + return rgbImage; + } else if (rgbImage.data.length === height * width * 3) { // Need to add alpha channel - const rgbaLength = (rgbImage.length / 3) * 4; + const rgbaLength = (rgbImage.data.length / 3) * 4; const rgbaArray = new Uint8ClampedArray(rgbaLength); - for (let i = 0; i < rgbImage.length / 3; ++i) { - rgbaArray[i * 4] = rgbImage[i * 3]!; - rgbaArray[i * 4 + 1] = rgbImage[i * 3 + 1]!; - rgbaArray[i * 4 + 2] = rgbImage[i * 3 + 2]!; + for (let i = 0; i < rgbImage.data.length / 3; ++i) { + rgbaArray[i * 4] = rgbImage.data[i * 3]!; + rgbaArray[i * 4 + 1] = rgbImage.data[i * 3 + 1]!; + rgbaArray[i * 4 + 2] = rgbImage.data[i * 3 + 2]!; rgbaArray[i * 4 + 3] = 255; } - return new ImageData(rgbaArray, width, height); + return { + ...rgbImage, + count: 4, + data: rgbaArray, + }; } else { throw new Error( - `Unexpected number of channels in raster data: ${rgbImage.length / (height * width)}`, + `Unexpected number of channels in raster data: ${rgbImage.data.length / (height * width)}`, ); } } @@ -137,24 +79,14 @@ export function parseColormap(cmap: Uint16Array): ImageData { } export async function fetchGeoTIFF( - input: GeoTIFF | string | ArrayBuffer | Blob | BaseClient, + input: GeoTIFF | string | URL | ArrayBuffer, ): Promise { - if (typeof input === "string") { - return fromUrl(input); + if (typeof input === "string" || input instanceof URL) { + return await GeoTIFF.fromUrl(input); } if (input instanceof ArrayBuffer) { - return fromArrayBuffer(input); - } - - if (input instanceof Blob) { - return fromBlob(input); - } - - // TODO: instanceof may fail here if multiple versions of geotiff.js are - // present - if (input instanceof BaseClient) { - return fromCustomClient(input); + return await GeoTIFF.fromArrayBuffer(input); } return input; @@ -164,18 +96,12 @@ export async function fetchGeoTIFF( * Calculate the WGS84 bounding box of a GeoTIFF image */ export function getGeographicBounds( - image: GeoTIFFImage, + geotiff: GeoTIFF, converter: Converter, ): { west: number; south: number; east: number; north: number } { - const projectedBbox = image.getBoundingBox() as [ - number, - number, - number, - number, - ]; + const [minX, minY, maxX, maxY] = geotiff.bbox; // Reproject all four corners to handle rotation/skew - const [minX, minY, maxX, maxY] = projectedBbox; const corners: [number, number][] = [ converter.forward([minX, minY]), // bottom-left converter.forward([maxX, minY]), // bottom-right @@ -195,20 +121,3 @@ export function getGeographicBounds( // Return bounds in MapLibre format: [[west, south], [east, north]] return { west, south, east, north }; } - -/** Parse the GDAL_NODATA TIFF tag into a number. */ -export function parseGDALNoData( - GDAL_NODATA: string | undefined, -): number | null { - if (!GDAL_NODATA) { - return null; - } - - // Remove trailing null character if present - const noDataString = - GDAL_NODATA?.[GDAL_NODATA?.length - 1] === "\x00" - ? GDAL_NODATA.slice(0, -1) - : GDAL_NODATA; - - return noDataString?.length > 0 ? parseFloat(noDataString) : null; -} diff --git a/packages/deck.gl-geotiff/src/geotiff/render-pipeline.ts b/packages/deck.gl-geotiff/src/geotiff/render-pipeline.ts index f3650ec..dca27f4 100644 --- a/packages/deck.gl-geotiff/src/geotiff/render-pipeline.ts +++ b/packages/deck.gl-geotiff/src/geotiff/render-pipeline.ts @@ -1,3 +1,4 @@ +import { Photometric, SampleFormat } from "@cogeotiff/core"; import type { RasterModule } from "@developmentseed/deck.gl-raster/gpu-modules"; import { CMYKToRGB, @@ -7,13 +8,11 @@ import { FilterNoDataVal, YCbCrToRGB, } from "@developmentseed/deck.gl-raster/gpu-modules"; +import type { GeoTIFF, Overview } from "@developmentseed/geotiff"; import type { Device, SamplerProps, Texture } from "@luma.gl/core"; -import type { GeoTIFFImage, TypedArrayWithDimensions } from "geotiff"; import type { COGLayerProps, GetTileDataOptions } from "../cog-layer"; -import { addAlphaChannel, parseColormap, parseGDALNoData } from "./geotiff"; +import { addAlphaChannel, parseColormap } from "./geotiff"; import { inferTextureFormat } from "./texture"; -import type { ImageFileDirectory } from "./types"; -import { PhotometricInterpretationT } from "./types"; export type TextureDataT = { height: number; @@ -41,19 +40,21 @@ type UnresolvedRasterModule = }; export function inferRenderPipeline( - // TODO: narrow type to only used fields - ifd: ImageFileDirectory, + geotiff: GeoTIFF, device: Device, ): { getTileData: COGLayerProps["getTileData"]; renderTile: COGLayerProps["renderTile"]; } { - const { SampleFormat } = ifd; + const { sampleFormat } = geotiff.cachedTags; + if (sampleFormat === null) { + throw new Error("SampleFormat tag is required to infer render pipeline"); + } - switch (SampleFormat[0]) { + switch (sampleFormat[0]) { // Unsigned integers - case 1: - return createUnormPipeline(ifd, device); + case SampleFormat.Uint: + return createUnormPipeline(geotiff, device); } throw new Error( @@ -65,20 +66,21 @@ export function inferRenderPipeline( * Create pipeline for visualizing unsigned-integer data. */ function createUnormPipeline( - ifd: ImageFileDirectory, + geotiff: GeoTIFF, device: Device, ): { getTileData: COGLayerProps["getTileData"]; renderTile: COGLayerProps["renderTile"]; } { + const tags = geotiff.cachedTags; const { - BitsPerSample, - ColorMap, - GDAL_NODATA, - PhotometricInterpretation, - SampleFormat, - SamplesPerPixel, - } = ifd; + bitsPerSample, + colorMap, + photometric, + sampleFormat, + samplesPerPixel, + nodata, + } = tags; const renderPipeline: UnresolvedRasterModule[] = [ { @@ -89,11 +91,9 @@ function createUnormPipeline( }, ]; - // Add NoData filtering if GDAL_NODATA is defined - const noDataVal = parseGDALNoData(GDAL_NODATA); - if (noDataVal !== null) { + if (nodata !== null) { // Since values are 0-1 for unorm textures, - const noDataScaled = noDataVal / 255.0; + const noDataScaled = nodata / 255.0; renderPipeline.push({ module: FilterNoDataVal, @@ -102,9 +102,9 @@ function createUnormPipeline( } const toRGBModule = photometricInterpretationToRGB( - PhotometricInterpretation, + photometric, device, - ColorMap, + colorMap, ); if (toRGBModule) { renderPipeline.push(toRGBModule); @@ -112,7 +112,7 @@ function createUnormPipeline( // For palette images, use nearest-neighbor sampling const samplerOptions: SamplerProps = - PhotometricInterpretation === PhotometricInterpretationT.Palette + photometric === Photometric.Palette ? { magFilter: "nearest", minFilter: "nearest", @@ -123,44 +123,44 @@ function createUnormPipeline( }; const getTileData: COGLayerProps["getTileData"] = async ( - image: GeoTIFFImage, + image: GeoTIFF | Overview, options: GetTileDataOptions, ) => { - const { device } = options; - const mergedOptions = { - ...options, - interleave: true, - }; + const { device, x, y } = options; + // TODO: pass down signal + const tile = await image.fetchTile(x, y); + let { array } = tile; - let data: TypedArrayWithDimensions | ImageData = (await image.readRasters( - mergedOptions, - )) as TypedArrayWithDimensions; - let numSamples = SamplesPerPixel; + let numSamples = samplesPerPixel; - if (SamplesPerPixel === 3) { + if (samplesPerPixel === 3) { // WebGL2 doesn't have an RGB-only texture format; it requires RGBA. - data = addAlphaChannel(data); + array = addAlphaChannel(array); numSamples = 4; } + if (array.layout === "band-separate") { + throw new Error("Band-separate images not yet implemented."); + } + const textureFormat = inferTextureFormat( // Add one sample for added alpha channel numSamples, - BitsPerSample, - SampleFormat, + bitsPerSample, + sampleFormat, ); const texture = device.createTexture({ - data, + data: array.data, format: textureFormat, - width: data.width, - height: data.height, + width: array.width, + height: array.height, sampler: samplerOptions, }); return { texture, - height: data.height, - width: data.width, + height: array.height, + width: array.width, }; }; const renderTile: COGLayerProps["renderTile"] = ( @@ -173,14 +173,14 @@ function createUnormPipeline( } function photometricInterpretationToRGB( - PhotometricInterpretation: number, + photometric: Photometric, device: Device, ColorMap?: Uint16Array, ): RasterModule | null { - switch (PhotometricInterpretation) { - case PhotometricInterpretationT.RGB: + switch (photometric) { + case Photometric.Rgb: return null; - case PhotometricInterpretationT.Palette: { + case Photometric.Palette: { if (!ColorMap) { throw new Error( "ColorMap is required for PhotometricInterpretation Palette", @@ -207,22 +207,21 @@ function photometricInterpretationToRGB( }; } - case PhotometricInterpretationT.CMYK: + // Not sure why cogeotiff calls this "Separated", but it means CMYK + case Photometric.Separated: return { module: CMYKToRGB, }; - case PhotometricInterpretationT.YCbCr: + case Photometric.Ycbcr: return { module: YCbCrToRGB, }; - case PhotometricInterpretationT.CIELab: + case Photometric.Cielab: return { module: cieLabToRGB, }; default: - throw new Error( - `Unsupported PhotometricInterpretation ${PhotometricInterpretation}`, - ); + throw new Error(`Unsupported PhotometricInterpretation ${photometric}`); } } diff --git a/packages/deck.gl-geotiff/src/geotiff/texture.ts b/packages/deck.gl-geotiff/src/geotiff/texture.ts index 2c881c3..f54d896 100644 --- a/packages/deck.gl-geotiff/src/geotiff/texture.ts +++ b/packages/deck.gl-geotiff/src/geotiff/texture.ts @@ -1,22 +1,21 @@ -import type { TextureFormat, TextureProps } from "@luma.gl/core"; -import type { GeoTIFFImage, TypedArray } from "geotiff"; -import type { ImageFileDirectory } from "./types"; +import { SampleFormat } from "@cogeotiff/core"; +import type { GeoTIFF } from "@developmentseed/geotiff"; +import type { TextureFormat, TextureProps, TypedArray } from "@luma.gl/core"; /** * Infers texture properties from a GeoTIFF image and its associated data. */ export function createTextureProps( - image: GeoTIFFImage, + geotiff: GeoTIFF, data: TypedArray, options: { width: number; height: number }, ): TextureProps { - const ifd = image.getFileDirectory() as ImageFileDirectory; - const samplesPerPixel = ifd.SamplesPerPixel; + const { samplesPerPixel, bitsPerSample, sampleFormat } = geotiff.cachedTags; const textureFormat = inferTextureFormat( samplesPerPixel, - ifd.BitsPerSample, - ifd.SampleFormat, + bitsPerSample, + sampleFormat, ); return { @@ -33,7 +32,7 @@ export function createTextureProps( export function inferTextureFormat( samplesPerPixel: number, bitsPerSample: Uint16Array, - sampleFormat: Uint16Array, + sampleFormat: SampleFormat[], ): TextureFormat { const channelCount = verifySamplesPerPixel(samplesPerPixel); const bitWidth = verifyIdenticalBitsPerSample(bitsPerSample); @@ -94,7 +93,7 @@ function verifyIdenticalBitsPerSample(bitsPerSample: Uint16Array): BitWidth { /** * Map the geotiff tag SampleFormat to known kinds of scalars */ -function inferScalarKind(sampleFormat: Uint16Array): ScalarKind { +function inferScalarKind(sampleFormat: SampleFormat[]): ScalarKind { // Only support identical SampleFormats for all samples const first = sampleFormat[0]!; @@ -107,11 +106,11 @@ function inferScalarKind(sampleFormat: Uint16Array): ScalarKind { } switch (first) { - case 1: + case SampleFormat.Uint: return "unorm"; - case 2: + case SampleFormat.Int: return "sint"; - case 3: + case SampleFormat.Float: return "float"; default: throw new Error(`Unsupported SampleFormat ${sampleFormat}`); diff --git a/packages/deck.gl-geotiff/src/geotiff/types.ts b/packages/deck.gl-geotiff/src/geotiff/types.ts deleted file mode 100644 index 995f19c..0000000 --- a/packages/deck.gl-geotiff/src/geotiff/types.ts +++ /dev/null @@ -1,35 +0,0 @@ -export enum PhotometricInterpretationT { - WhiteIsZero = 0, - BlackIsZero = 1, - RGB = 2, - Palette = 3, - TransparencyMask = 4, - CMYK = 5, - YCbCr = 6, - CIELab = 8, - ICCLab = 9, -} - -export enum PlanarConfigurationT { - Chunky = 1, - Planar = 2, -} - -/** Improved typing for IFD. */ -export type ImageFileDirectory = { - BitsPerSample: Uint16Array; - ColorMap?: Uint16Array; - Compression: number; - /** GDAL NoData value as string. - * - */ - GDAL_NODATA?: string; - ImageLength: number; - ImageWidth: number; - PhotometricInterpretation: PhotometricInterpretationT; - /** Strip or tiled */ - PlanarConfiguration: PlanarConfigurationT; - SampleFormat: Uint16Array; - /** Number of bands */ - SamplesPerPixel: number; -}; diff --git a/packages/deck.gl-geotiff/src/index.ts b/packages/deck.gl-geotiff/src/index.ts index fe8e598..bd6e396 100644 --- a/packages/deck.gl-geotiff/src/index.ts +++ b/packages/deck.gl-geotiff/src/index.ts @@ -1,14 +1,10 @@ export type { COGLayerProps } from "./cog-layer.js"; export { COGLayer } from "./cog-layer.js"; -export { parseCOGTileMatrixSet } from "./cog-tile-matrix-set.js"; -export { loadRgbImage, parseColormap } from "./geotiff/geotiff.js"; +export { parseColormap } from "./geotiff/geotiff.js"; export * as texture from "./geotiff/texture.js"; -export type { GeoTIFFLayerProps } from "./geotiff-layer.js"; -export { GeoTIFFLayer } from "./geotiff-layer.js"; -export { - extractGeotiffReprojectors, - fromGeoTransform, -} from "./geotiff-reprojection.js"; +// Don't export GeoTIFF Layer for now; nudge people towards COGLayer +// export type { GeoTIFFLayerProps } from "./geotiff-layer.js"; +// export { GeoTIFFLayer } from "./geotiff-layer.js"; export type { MosaicLayerProps } from "./mosaic-layer/mosaic-layer.js"; export { MosaicLayer } from "./mosaic-layer/mosaic-layer.js"; export { @@ -16,7 +12,3 @@ export { MosaicTileset2D, } from "./mosaic-layer/mosaic-tileset-2d"; export * as proj from "./proj.js"; - -import { __TEST_EXPORTS as cogTileMatrixSetTestExports } from "./cog-tile-matrix-set.js"; - -export const __TEST_EXPORTS = { ...cogTileMatrixSetTestExports }; diff --git a/packages/deck.gl-geotiff/src/proj.ts b/packages/deck.gl-geotiff/src/proj.ts index 95f00a5..077ad3b 100644 --- a/packages/deck.gl-geotiff/src/proj.ts +++ b/packages/deck.gl-geotiff/src/proj.ts @@ -1,82 +1,37 @@ -import type { ProjectionDefinition } from "proj4"; -import proj4 from "proj4"; -import type { PROJJSONDefinition } from "proj4/dist/lib/core"; - -export type SupportedCrsUnit = - | "m" - | "metre" - | "meter" - | "meters" - | "foot" - | "US survey foot" - | "degree"; - -export interface ProjectionInfo { - /** Proj4-compatible projection definition (PROJJSON or proj4 string) */ - def: string | PROJJSONDefinition; - /** A parsed projection definition */ - parsed: ProjectionDefinition; - /** Units of the coordinate system */ - coordinatesUnits: SupportedCrsUnit; -} - -export type GeoKeysParser = ( - geoKeys: Record, -) => Promise; +import type { ProjJson } from "@developmentseed/geotiff"; +import type { ProjectionDefinition } from "wkt-parser"; +import wktParser from "wkt-parser"; /** - * Get the Projection of a GeoTIFF - * - * The first `image` must be passed in, as only the top-level IFD contains geo - * keys. + * A global registry holding parsed projection definitions. */ -export async function epsgIoGeoKeyParser( - geoKeys: Record, -): Promise { - const projectionCode: number | null = - geoKeys.ProjectedCSTypeGeoKey || geoKeys.GeographicTypeGeoKey || null; +const PROJECTION_REGISTRY = new Map(); - const sourceProjection = await getProjjson(projectionCode); +export type EpsgResolver = (epsg: number) => Promise; - if (!sourceProjection) { - return null; +export async function epsgResolver(epsg: number) { + const key = `EPSG:${epsg}`; + // TODO: override global proj4 EPSG:4326 definition because it doesn't include + // semi major axis + const cachedProj = PROJECTION_REGISTRY.get(key); + if (cachedProj !== undefined) { + return cachedProj; } - const parsed = parseCrs(sourceProjection); - return { - def: sourceProjection, - parsed, - coordinatesUnits: parsed.units as SupportedCrsUnit, - }; + const projjson = await getProjjson(epsg); + const proj = wktParser(projjson); + PROJECTION_REGISTRY.set(key, proj); + + return proj; } /** Query epsg.io for the PROJJSON corresponding to the given EPSG code. */ -async function getProjjson(projectionCode: number | null) { - if (projectionCode === null) { - return null; +async function getProjjson(epsg: number): Promise { + const url = `https://epsg.io/${epsg}.json`; + const resp = await fetch(url); + if (!resp.ok) { + throw new Error(`Failed to fetch PROJJSON from ${url}`); } - 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; -} - -/** - * Parse a proj4-accepted input into a ProjectionDefinition - */ -export function parseCrs( - crs: string | PROJJSONDefinition, -): ProjectionDefinition { - // If you pass proj4.defs a projjson, it doesn't parse it; it just returns the - // input. - // - // Instead, you need to assign it to an alias and then retrieve it. - - const key = "__deck.gl-geotiff-internal__"; - proj4.defs(key, crs); - return proj4.defs(key); + return await resp.json(); } diff --git a/packages/deck.gl-geotiff/src/wkt-parser.d.ts b/packages/deck.gl-geotiff/src/wkt-parser.d.ts new file mode 100644 index 0000000..8573bd1 --- /dev/null +++ b/packages/deck.gl-geotiff/src/wkt-parser.d.ts @@ -0,0 +1,62 @@ +declare module "wkt-parser" { + import type { ProjJson } from "./crs.js"; + + export interface DatumDefinition { + /** The type of datum. */ + datum_type: number; + /** Semi-major axis of the ellipsoid. */ + a: number; + /** Semi-minor axis of the ellipsoid. */ + b: number; + /** Eccentricity squared of the ellipsoid. */ + es: number; + /** Second eccentricity squared of the ellipsoid. */ + ep2: number; + } + + export interface ProjectionDefinition { + title: string; + projName?: string; + ellps?: string; + datum?: DatumDefinition; + datumName?: string; + rf?: number; + lat0?: number; + lat1?: number; + lat2?: number; + lat_ts?: number; + long0?: number; + long1?: number; + long2?: number; + alpha?: number; + longc?: number; + x0?: number; + y0?: number; + k0?: number; + a?: number; + b?: number; + R_A?: true; + zone?: number; + utmSouth?: true; + datum_params?: string | number[]; + to_meter?: number; + units?: string; + from_greenwich?: number; + datumCode?: string; + nadgrids?: string; + axis?: string; + sphere?: boolean; + rectified_grid_angle?: number; + approx?: boolean; + over?: boolean; + projStr?: string; + } + + /** + * Parse a WKT string or PROJJSON object into a proj4-compatible projection + * definition. + */ + export default function wktParser( + input: string | ProjJson, + ): ProjectionDefinition; +} diff --git a/packages/deck.gl-geotiff/tests/cog-tms.test.ts b/packages/deck.gl-geotiff/tests/cog-tms.test.ts deleted file mode 100644 index ee4f9c0..0000000 --- a/packages/deck.gl-geotiff/tests/cog-tms.test.ts +++ /dev/null @@ -1,99 +0,0 @@ -import { fromUrl } from "geotiff"; -import { describe, expect, it } from "vitest"; -import { parseCOGTileMatrixSet, __TEST_EXPORTS as TestExports } from "../src"; -import { epsgIoGeoKeyParser } from "../src/proj"; - -describe("create TileMatrixSet from COG", () => { - it("creates TMS", async () => { - const url = - "https://ds-wheels.s3.us-east-1.amazonaws.com/m_4007307_sw_18_060_20220803.tif"; - - const tiff = await fromUrl(url); - const tms = await parseCOGTileMatrixSet(tiff, epsgIoGeoKeyParser); - - const expectedTileMatrices = [ - { - id: "0", - scaleDenominator: 68684.95742667928, - cellSize: 19.231788079470196, - pointOfOrigin: [647118, 4533600], - tileWidth: 512, - tileHeight: 512, - matrixWidth: 1, - matrixHeight: 1, - geotransform: [ - 19.231788079470196, 0, 647118, 0, -19.231788079470196, 4533600, - ], - }, - { - id: "1", - scaleDenominator: 34285.71428571429, - cellSize: 9.6, - pointOfOrigin: [647118, 4533600], - tileWidth: 512, - tileHeight: 512, - matrixWidth: 2, - matrixHeight: 2, - geotransform: [9.6, 0, 647118, 0, -9.6, 4533600], - }, - { - id: "2", - scaleDenominator: 17142.857142857145, - cellSize: 4.8, - pointOfOrigin: [647118, 4533600], - tileWidth: 512, - tileHeight: 512, - matrixWidth: 3, - matrixHeight: 4, - geotransform: [4.8, 0, 647118, 0, -4.8, 4533600], - }, - { - id: "3", - scaleDenominator: 8571.428571428572, - cellSize: 2.4, - pointOfOrigin: [647118, 4533600], - tileWidth: 512, - tileHeight: 512, - matrixWidth: 5, - matrixHeight: 7, - geotransform: [2.4, 0, 647118, 0, -2.4, 4533600], - }, - { - id: "4", - scaleDenominator: 4285.714285714286, - cellSize: 1.2, - pointOfOrigin: [647118, 4533600], - tileWidth: 512, - tileHeight: 512, - matrixWidth: 10, - matrixHeight: 13, - geotransform: [1.2, 0, 647118, 0, -1.2, 4533600], - }, - { - id: "5", - scaleDenominator: 2142.857142857143, - cellSize: 0.6, - pointOfOrigin: [647118, 4533600], - tileWidth: 512, - tileHeight: 512, - matrixWidth: 19, - matrixHeight: 25, - geotransform: [0.6, 0, 647118, 0, -0.6, 4533600], - }, - ]; - - expect(tms.tileMatrices).toStrictEqual(expectedTileMatrices); - }); -}); - -describe("metersPerUnit", () => { - it("handles lowercase us survey foot", () => { - // @ts-expect-error testing case insensitivity with standard casing - expect(TestExports.metersPerUnit({}, "us survey foot")).toBe(1200 / 3937); - }); - - it("handles mixed case US Survey Foot", () => { - // @ts-expect-error testing case insensitivity with non-standard casing - expect(TestExports.metersPerUnit({}, "US Survey Foot")).toBe(1200 / 3937); - }); -}); diff --git a/packages/deck.gl-geotiff/tests/render-pipeline.test.ts b/packages/deck.gl-geotiff/tests/render-pipeline.test.ts index 11fab5f..20058b3 100644 --- a/packages/deck.gl-geotiff/tests/render-pipeline.test.ts +++ b/packages/deck.gl-geotiff/tests/render-pipeline.test.ts @@ -1,8 +1,9 @@ +import { Photometric, SampleFormat } from "@cogeotiff/core"; import type { RasterModule } from "@developmentseed/deck.gl-raster"; -import { globals } from "geotiff"; +import type { GeoTIFF } from "@developmentseed/geotiff"; import { describe, expect, it } from "vitest"; +import type { CachedTags } from "../../geotiff/dist/ifd"; import { inferRenderPipeline } from "../src/geotiff/render-pipeline"; -import type { ImageFileDirectory } from "../src/geotiff/types"; const MOCK_DEVICE = { createTexture: (x: any) => x, @@ -11,21 +12,20 @@ const MOCK_RENDER_TILE_DATA = { texture: {}, }; -// import {} from "@" type RelevantImageFileDirectory = Pick< - ImageFileDirectory, - | "BitsPerSample" - | "ColorMap" - | "GDAL_NODATA" - | "PhotometricInterpretation" - | "SampleFormat" - | "SamplesPerPixel" + CachedTags, + | "bitsPerSample" + | "colorMap" + | "nodata" + | "photometric" + | "sampleFormat" + | "samplesPerPixel" >; describe("land cover, single-band uint8", () => { const ifd: RelevantImageFileDirectory = { - BitsPerSample: new Uint16Array([8]), - ColorMap: new Uint16Array([ + bitsPerSample: new Uint16Array([8]), + colorMap: new Uint16Array([ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 17990, 53713, 0, 0, 0, 0, 0, 0, 0, 0, 57054, 55769, 60395, 43947, 0, 0, 0, 0, 0, 0, 46003, 0, 0, 0, 0, 0, 0, 0, 0, 0, 26728, 7196, 46517, 0, 0, 0, 0, 0, 0, 0, 0, 52428, 0, 0, 0, 0, 0, 0, @@ -61,14 +61,17 @@ describe("land cover, single-band uint8", () => { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ]), - GDAL_NODATA: "250\u0000", - PhotometricInterpretation: globals.photometricInterpretations.Palette, - SampleFormat: new Uint16Array([1]), - SamplesPerPixel: 1, + nodata: 250, + photometric: Photometric.Palette, + sampleFormat: [SampleFormat.Uint], + samplesPerPixel: 1, }; + const geotiff = { + cachedTags: ifd, + } as unknown as GeoTIFF; const { getTileData: _, renderTile } = inferRenderPipeline( - ifd as ImageFileDirectory, + geotiff, MOCK_DEVICE as any, ); const renderPipeline = renderTile( diff --git a/packages/deck.gl-geotiff/tsconfig.build.json b/packages/deck.gl-geotiff/tsconfig.build.json index 86953b5..546e11e 100644 --- a/packages/deck.gl-geotiff/tsconfig.build.json +++ b/packages/deck.gl-geotiff/tsconfig.build.json @@ -8,7 +8,9 @@ "include": ["src/**/*"], "exclude": ["node_modules", "dist", "tests"], "references": [ - { "path": "../raster-reproject/tsconfig.build.json" }, - { "path": "../deck.gl-raster/tsconfig.build.json" } + { "path": "../deck.gl-raster/tsconfig.build.json" }, + { "path": "../geotiff/tsconfig.build.json" }, + { "path": "../morecantile/tsconfig.build.json" }, + { "path": "../raster-reproject/tsconfig.build.json" } ] } diff --git a/packages/deck.gl-raster/package.json b/packages/deck.gl-raster/package.json index a23fdab..7204f48 100644 --- a/packages/deck.gl-raster/package.json +++ b/packages/deck.gl-raster/package.json @@ -60,6 +60,8 @@ "@luma.gl/shadertools": "^9.2.7" }, "dependencies": { + "@developmentseed/affine": "workspace:^", + "@developmentseed/morecantile": "workspace:^", "@developmentseed/raster-reproject": "workspace:^", "@math.gl/core": "^4.1.0", "@math.gl/culling": "^4.1.0", diff --git a/packages/deck.gl-raster/src/index.ts b/packages/deck.gl-raster/src/index.ts index 1189835..f8127f4 100644 --- a/packages/deck.gl-raster/src/index.ts +++ b/packages/deck.gl-raster/src/index.ts @@ -1,12 +1,7 @@ export type { RasterModule } from "./gpu-modules/types.js"; export type { RasterLayerProps } from "./raster-layer.js"; export { RasterLayer } from "./raster-layer.js"; -export { RasterTileset2D } from "./raster-tileset/index.js"; -export type { - TileMatrix, - TileMatrixSet, - TileMatrixSetBoundingBox, -} from "./raster-tileset/types.js"; +export { TileMatrixSetTileset } from "./raster-tileset/index.js"; import { __TEST_EXPORTS as traversalTestExports } from "./raster-tileset/raster-tile-traversal.js"; diff --git a/packages/deck.gl-raster/src/raster-tileset/index.ts b/packages/deck.gl-raster/src/raster-tileset/index.ts index d19311b..f96ab6b 100644 --- a/packages/deck.gl-raster/src/raster-tileset/index.ts +++ b/packages/deck.gl-raster/src/raster-tileset/index.ts @@ -1,2 +1 @@ -export { RasterTileset2D } from "./raster-tileset-2d.js"; -export type { TileMatrix, TileMatrixSet } from "./types.js"; +export { TileMatrixSetTileset } from "./raster-tileset-2d.js"; 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 9251262..fbcbe02 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 @@ -17,6 +17,8 @@ import type { Viewport } from "@deck.gl/core"; import { _GlobeViewport, assert } from "@deck.gl/core"; +import type { TileMatrix, TileMatrixSet } from "@developmentseed/morecantile"; +import { xy_bounds } from "@developmentseed/morecantile"; import type { OrientedBoundingBox } from "@math.gl/culling"; import { CullingVolume, @@ -27,9 +29,9 @@ import { lngLatToWorld, worldToLngLat } from "@math.gl/web-mercator"; import type { Bounds, + CornerBounds, + ProjectionFunction, TileIndex, - TileMatrix, - TileMatrixSet, ZRange, } from "./types.js"; @@ -138,11 +140,22 @@ export class RasterTileNode { /** A cache of the children of this node. */ private _children?: RasterTileNode[] | null; - constructor(x: number, y: number, z: number, metadata: TileMatrixSet) { + private projectTo3857: ProjectionFunction; + + constructor( + x: number, + y: number, + z: number, + { + metadata, + projectTo3857, + }: { metadata: TileMatrixSet; projectTo3857: ProjectionFunction }, + ) { this.x = x; this.y = y; this.z = z; this.metadata = metadata; + this.projectTo3857 = projectTo3857; } /** Get overview info for this tile's z level */ @@ -172,12 +185,9 @@ export class RasterTileNode { const childMatrix = this.metadata.tileMatrices[childZ]!; // Compute this tile's bounds in TMS' CRS - const parentBounds = computeProjectedTileBounds({ + const parentBounds = computeProjectedTileBounds(parentMatrix, { x: this.x, y: this.y, - transform: parentMatrix.geotransform, - tileWidth: parentMatrix.tileWidth, - tileHeight: parentMatrix.tileHeight, }); // Find overlapping child index range @@ -188,9 +198,15 @@ export class RasterTileNode { const children: RasterTileNode[] = []; + const { metadata, projectTo3857 } = this; for (let y = minRow; y <= maxRow; y++) { for (let x = minCol; x <= maxCol; x++) { - children.push(new RasterTileNode(x, y, childZ, this.metadata)); + children.push( + new RasterTileNode(x, y, childZ, { + metadata, + projectTo3857, + }), + ); } } @@ -364,12 +380,9 @@ export class RasterTileNode { const [minX, minY, maxX, maxY] = bounds; const [tileMinX, tileMinY, tileMaxX, tileMaxY] = commonSpaceBounds; - // console.log("bounds:", bounds); - // console.log("tile bounds:", commonSpaceBounds); - const inside = tileMinX < maxX && tileMaxX > minX && tileMinY < maxY && tileMaxY > minY; - // console.log("insideBounds", inside); + return inside; } @@ -414,21 +427,17 @@ export class RasterTileNode { commonSpaceBounds: Bounds; } { const tileMatrix = this.tileMatrix; - const { tileWidth, tileHeight, geotransform } = tileMatrix; const [minZ, maxZ] = zRange; - const tileCrsBounds = computeProjectedTileBounds({ + const tileCrsBounds = computeProjectedTileBounds(tileMatrix, { x: this.x, y: this.y, - transform: geotransform, - tileWidth, - tileHeight, }); const refPointsEPSG3857 = sampleReferencePointsInEPSG3857( REF_POINTS_9, tileCrsBounds, - this.metadata.projectTo3857, + this.projectTo3857, ); const commonSpacePositions = refPointsEPSG3857.map((xy) => @@ -476,60 +485,22 @@ export class RasterTileNode { * * @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; - - // 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. +function computeProjectedTileBounds( + tileMatrix: TileMatrix, + { + x, + y, + }: { + x: number; + y: number; + }, +): [number, number, number, number] { + const bounds = xy_bounds(tileMatrix, { x, y }); return [ - Math.min(minX, maxX), - Math.min(minY, maxY), - Math.max(minX, maxX), - Math.max(minY, maxY), + bounds.lowerLeft[0], + bounds.lowerLeft[1], + bounds.upperRight[0], + bounds.upperRight[1], ]; } @@ -546,7 +517,7 @@ function computeProjectedTileBounds({ function sampleReferencePointsInEPSG3857( refPoints: [number, number][], tileBounds: [number, number, number, number], - projectTo3857: (xy: [number, number]) => [number, number], + projectTo3857: ProjectionFunction, ): [number, number][] { const [minX, minY, maxX, maxY] = tileBounds; const refPointPositions: [number, number][] = []; @@ -556,7 +527,7 @@ function sampleReferencePointsInEPSG3857( const geoY = minY + relY * (maxY - minY); // Reproject to Web Mercator (EPSG 3857) - const projected = projectTo3857([geoX, geoY]); + const projected = projectTo3857(geoX, geoY); refPointPositions.push(projected); } @@ -675,9 +646,11 @@ export function getTileIndices( viewport: Viewport; maxZ: number; zRange: ZRange | null; + projectTo3857: ProjectionFunction; + wgs84Bounds: CornerBounds; }, ): TileIndex[] { - const { viewport, maxZ, zRange } = opts; + const { viewport, maxZ, zRange, wgs84Bounds } = opts; // Only define `project` function for Globe viewports, same as upstream const project: ((xyz: number[]) => number[]) | null = @@ -723,7 +696,7 @@ export function getTileIndices( // minZ to 0 const minZ = 0; - const { lowerLeft, upperRight } = metadata.wgsBounds; + const { lowerLeft, upperRight } = wgs84Bounds; const [minLng, minLat] = lowerLeft; const [maxLng, maxLat] = upperRight; const bottomLeft = lngLatToWorld([minLng, minLat]); @@ -744,7 +717,12 @@ export function getTileIndices( const roots: RasterTileNode[] = []; for (let y = 0; y < rootMatrix.matrixHeight; y++) { for (let x = 0; x < rootMatrix.matrixWidth; x++) { - roots.push(new RasterTileNode(x, y, 0, metadata)); + roots.push( + new RasterTileNode(x, y, 0, { + metadata, + projectTo3857: opts.projectTo3857, + }), + ); } } @@ -759,9 +737,6 @@ export function getTileIndices( bounds, }; - // console.log("traversalParams", traversalParams); - // console.log("roots", roots); - for (const root of roots) { root.update(traversalParams); } @@ -800,19 +775,6 @@ function getMetersPerPixelAtBoundingVolume( return getMetersPerPixel(lat, zoom); } -// function getScreenMetersPerPixel(viewport: Viewport, center: Vector3): number { -// const lng -// const p0 = viewport.projectPosition(center); -// const p1 = viewport.projectPosition([ -// centerArray[0]! + 1, -// centerArray[1]!, -// centerArray[2]!, -// ]); - -// const pixelsPerMeter = Math.hypot(p1[0] - p0[0], p1[1] - p0[1]); -// return 1 / pixelsPerMeter; -// } - /** * Exports only for use in testing */ 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 c10e0b3..8d3a3b6 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 @@ /** - * RasterTileset2D - Improved Implementation with Frustum Culling + * TileMatrixSetTileset - 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. @@ -8,20 +8,55 @@ import type { Viewport } from "@deck.gl/core"; import type { _Tileset2DProps as Tileset2DProps } from "@deck.gl/geo-layers"; import { _Tileset2D as Tileset2D } from "@deck.gl/geo-layers"; +import * as affine from "@developmentseed/affine"; +import type { BoundingBox, TileMatrixSet } from "@developmentseed/morecantile"; +import { tileTransform } from "@developmentseed/morecantile"; import type { Matrix4 } from "@math.gl/core"; import { getTileIndices } from "./raster-tile-traversal"; -import type { Bounds, TileIndex, TileMatrixSet, ZRange } from "./types"; +import type { + Bounds, + CornerBounds, + Point, + ProjectionFunction, + TileIndex, + ZRange, +} from "./types"; /** - * RasterTileset2D with proper frustum culling + * A generic tileset implementation organized according to the OGC + * [TileMatrixSet](https://docs.ogc.org/is/17-083r4/17-083r4.html) + * specification. */ -export class RasterTileset2D extends Tileset2D { - private metadata: TileMatrixSet; - - constructor(metadata: TileMatrixSet, opts: Tileset2DProps) { +export class TileMatrixSetTileset extends Tileset2D { + private tms: TileMatrixSet; + private wgs84Bounds: CornerBounds; + private projectTo3857: ProjectionFunction; + + constructor( + opts: Tileset2DProps, + tms: TileMatrixSet, + { + projectTo4326, + projectTo3857, + }: { + projectTo4326: ProjectionFunction; + projectTo3857: ProjectionFunction; + }, + ) { super(opts); - this.metadata = metadata; + this.tms = tms; + this.projectTo3857 = projectTo3857; + + if (!tms.boundingBox) { + throw new Error( + "Bounding Box inference not yet implemented; should be provided on TileMatrixSet", + ); + } + + this.wgs84Bounds = projectBoundsToWgs84(tms.boundingBox, projectTo4326, { + densifyPts: 10, + }); } /** @@ -38,20 +73,21 @@ export class RasterTileset2D extends Tileset2D { modelMatrix?: Matrix4; modelMatrixInverse?: Matrix4; }): TileIndex[] { - // console.log("Called getTileIndices", opts); - const maxAvailableZ = this.metadata.tileMatrices.length - 1; + const maxAvailableZ = this.tms.tileMatrices.length - 1; const maxZ = typeof opts.maxZoom === "number" ? Math.min(opts.maxZoom, maxAvailableZ) : maxAvailableZ; - const tileIndices = getTileIndices(this.metadata, { + const tileIndices = getTileIndices(this.tms, { viewport: opts.viewport, maxZ, zRange: opts.zRange ?? null, + wgs84Bounds: this.wgs84Bounds, + projectTo3857: this.projectTo3857, }); - // console.log("Visible tile indices:", tileIndices.length); + return tileIndices; } @@ -65,8 +101,8 @@ 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.tms.tileMatrices[index.z]!; + const parentOverview = this.tms.tileMatrices[index.z - 1]!; const decimation = currentOverview.cellSize / parentOverview.cellSize; @@ -83,15 +119,10 @@ export class RasterTileset2D extends Tileset2D { override getTileMetadata(index: TileIndex): Record { const { x, y, z } = index; - const { tileMatrices } = this.metadata; + const { tileMatrices } = this.tms; 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] = geotransform; + const { tileHeight, tileWidth } = tileMatrix; + const tileAffine = tileTransform(tileMatrix, { col: x, row: y }); // Calculate pixel coordinates for this tile's extent const pixelMinCol = x * tileWidth; @@ -100,22 +131,10 @@ export class RasterTileset2D extends Tileset2D { const pixelMaxRow = (y + 1) * tileHeight; // Calculate the four corners of the tile in geographic coordinates - const topLeft: [number, number] = [ - a * pixelMinCol + b * pixelMinRow + c, - d * pixelMinCol + e * pixelMinRow + f, - ]; - const topRight: [number, number] = [ - a * pixelMaxCol + b * pixelMinRow + c, - d * pixelMaxCol + e * pixelMinRow + f, - ]; - const bottomLeft: [number, number] = [ - a * pixelMinCol + b * pixelMaxRow + c, - d * pixelMinCol + e * pixelMaxRow + f, - ]; - const bottomRight: [number, number] = [ - a * pixelMaxCol + b * pixelMaxRow + c, - d * pixelMaxCol + e * pixelMaxRow + f, - ]; + const topLeft = affine.apply(tileAffine, pixelMinCol, pixelMinRow); + const topRight = affine.apply(tileAffine, pixelMaxCol, pixelMinRow); + const bottomLeft = affine.apply(tileAffine, pixelMinCol, pixelMaxRow); + const bottomRight = affine.apply(tileAffine, pixelMaxCol, pixelMaxRow); // Return the projected bounds as four corners // This preserves rotation/skew information @@ -143,3 +162,54 @@ export class RasterTileset2D extends Tileset2D { }; } } + +function projectBoundsToWgs84( + bounds: BoundingBox, + projectTo4326: ProjectionFunction, + { densifyPts }: { densifyPts: number }, +): CornerBounds { + const { lowerLeft, upperRight } = bounds; + + // Four corners of the bounding box + const corners: Point[] = [ + lowerLeft, + [upperRight[0], lowerLeft[1]], + upperRight, + [lowerLeft[0], upperRight[1]], + ]; + + // Densify edges: interpolate densifyPts points along each edge + const points: Point[] = []; + for (let i = 0; i < corners.length; i++) { + const from = corners[i]!; + const to = corners[(i + 1) % corners.length]!; + // Include the start corner and all intermediate points (end corner + // will be included as the start of the next edge) + for (let j = 0; j <= densifyPts; j++) { + const t = j / (densifyPts + 1); + points.push([ + from[0] + (to[0] - from[0]) * t, + from[1] + (to[1] - from[1]) * t, + ]); + } + } + + // Reproject all points to WGS84 and compute the bounding box + let wgsMinX = Infinity; + let wgsMinY = Infinity; + let wgsMaxX = -Infinity; + let wgsMaxY = -Infinity; + + for (const [x, y] of points) { + const [lon, lat] = projectTo4326(x, y); + if (lon < wgsMinX) wgsMinX = lon; + if (lat < wgsMinY) wgsMinY = lat; + if (lon > wgsMaxX) wgsMaxX = lon; + if (lat > wgsMaxY) wgsMaxY = lat; + } + + return { + lowerLeft: [wgsMinX, wgsMinY], + upperRight: [wgsMaxX, wgsMaxY], + }; +} diff --git a/packages/deck.gl-raster/src/raster-tileset/types.ts b/packages/deck.gl-raster/src/raster-tileset/types.ts index c44c30b..597d470 100644 --- a/packages/deck.gl-raster/src/raster-tileset/types.ts +++ b/packages/deck.gl-raster/src/raster-tileset/types.ts @@ -32,175 +32,21 @@ export type Point = [number, number]; type CRS = any; +export type ProjectionFunction = (x: number, y: number) => Point; + /** - * Minimum bounding rectangle surrounding a 2D resource in the CRS indicated elsewhere + * Bounding box defined by two named corners */ -export interface TileMatrixSetBoundingBox { +export type CornerBounds = { lowerLeft: Point; upperRight: Point; - crs?: CRS; -} - -/** - * Represents a single resolution level in a raster tileset. - * - * COGs contain multiple resolution levels (overviews) for efficient - * visualization at different zoom levels. - * - * IMPORTANT: Overviews are ordered according to TileMatrixSet specification: - * - Index 0: Coarsest resolution (most zoomed out) - * - Index N: Finest resolution (most zoomed in) - * - * This matches the natural ordering where z increases with detail. - */ -export type TileMatrix = { - id: string; - - /** - * Scale denominator of this tile matrix. - * - * Defined as cellSize (meters per pixel) * meters per unit / 0.00028 - */ - scaleDenominator: number; - - /** - * Cell size of this tile matrix. - * - * CRS units per pixel (not necessarily meters). - */ - cellSize: number; - - // /** - // * Overview index in the TileMatrixSet ordering. - // * - Index 0: Coarsest resolution (most zoomed out) - // * - Higher indices: Progressively finer resolution - // * - // * This is the index in the COGMetadata.overviews array and represents - // * the natural ordering from coarse to fine. - // * - // * Note: This is different from GeoTIFF's internal level numbering where - // * level 0 is the full resolution image. - // * - // * @example - // * // For a COG with 4 resolutions: - // * index: 0 // Coarsest: 1250x1000 pixels (8x downsampled) - // * index: 1 // Medium: 2500x2000 pixels (4x downsampled) - // * index: 2 // Fine: 5000x4000 pixels (2x downsampled) - // * index: 3 // Finest: 10000x8000 pixels (full resolution) - // */ - // z: number; - - pointOfOrigin: Point; - - /** - * Width of each tile of this tile matrix in pixels. - */ - tileWidth: number; - /** - * Height of each tile of this tile matrix in pixels. - */ - tileHeight: number; - - /** - * Number of tiles in the X (horizontal) direction at this overview level. - * - * Calculated as: Math.ceil(width / tileWidth) - * - * @example - * // If tileWidth = 512: - * tilesX: 3 // z=0: ceil(1250 / 512) - * tilesX: 5 // z=1: ceil(2500 / 512) - * tilesX: 10 // z=2: ceil(5000 / 512) - * tilesX: 20 // z=3: ceil(10000 / 512) - */ - matrixWidth: number; - - /** - * Number of tiles in the Y (vertical) direction at this overview level. - * - * Calculated as: Math.ceil(height / tileHeight) - * - * @example - * // If tileHeight = 512: - * tilesY: 2 // z=0: ceil(1000 / 512) - * tilesY: 4 // z=1: ceil(2000 / 512) - * tilesY: 8 // z=2: ceil(4000 / 512) - * tilesY: 16 // z=3: ceil(8000 / 512) - */ - matrixHeight: number; - - // /** - // * Downsampling scale factor relative to full resolution (finest level). - // * - // * Indicates how much this overview is downsampled compared to the finest resolution. - // * - Scale factor of 1: Full resolution (finest level) - // * - Scale factor of 2: Half resolution - // * - Scale factor of 4: Quarter resolution - // * - Scale factor of 8: Eighth resolution (coarsest in this example) - // * - // * Common pattern: Each overview is 2x downsampled from the next finer level, - // * so scale factors are powers of 2: 8, 4, 2, 1 (from coarsest to finest) - // * - // * @example - // * scaleFactor: 8 // z=0: 1250x1000 (8x downsampled from finest) - // * scaleFactor: 4 // z=1: 2500x2000 (4x downsampled) - // * scaleFactor: 2 // z=2: 5000x4000 (2x downsampled) - // * scaleFactor: 1 // z=3: 10000x8000 (full resolution) - // */ - // scaleFactor: number; - - /** - * Affine geotransform for this overview level. - * - * Uses Python `affine` package ordering (NOT GDAL ordering): - * [a, b, c, d, e, f] where: - * - x_geo = a * col + b * row + c - * - y_geo = d * col + e * row + f - * - * Parameters: - * - a: pixel width (x resolution) - * - b: row rotation (typically 0) - * - c: x-coordinate of upper-left corner of the upper-left pixel - * - d: column rotation (typically 0) - * - e: pixel height (y resolution, typically negative) - * - f: y-coordinate of upper-left corner of the upper-left pixel - * - * @example - * // For a UTM image with 30m pixels: - * [30, 0, 440720, 0, -30, 3751320] - * // x_geo = 30 * col + 440720 - * // y_geo = -30 * row + 3751320 - */ - geotransform: [number, number, number, number, number, number]; }; /** - * COG Metadata extracted from GeoTIFF + * Minimum bounding rectangle surrounding a 2D resource in the CRS indicated elsewhere */ -export type TileMatrixSet = { - /** - * Title of this tile matrix set, normally used for display to a human - */ - title?: string; - /** - * Brief narrative description of this tile matrix set, normally available for display to a human - */ - description?: string; - crs: CRS; - - boundingBox?: TileMatrixSetBoundingBox; - /** - * Describes scale levels and its tile matrices - */ - tileMatrices: TileMatrix[]; - - /** - * Bounding box of this TMS in WGS84 lon/lat. - */ - wgsBounds: TileMatrixSetBoundingBox; - - projectToWgs84: (point: [number, number]) => [number, number]; - projectTo3857: (point: [number, number]) => [number, number]; +export type TileMatrixSetBoundingBox = CornerBounds & { + crs?: CRS; }; /** diff --git a/packages/deck.gl-raster/tsconfig.build.json b/packages/deck.gl-raster/tsconfig.build.json index 9ab31b2..0e89e61 100644 --- a/packages/deck.gl-raster/tsconfig.build.json +++ b/packages/deck.gl-raster/tsconfig.build.json @@ -7,5 +7,8 @@ }, "include": ["src/**/*"], "exclude": ["node_modules", "dist", "tests"], - "references": [{ "path": "../raster-reproject/tsconfig.build.json" }] + "references": [ + { "path": "../morecantile/tsconfig.build.json" }, + { "path": "../raster-reproject/tsconfig.build.json" } + ] } diff --git a/packages/geotiff/README.md b/packages/geotiff/README.md index 10b1b06..cdcd0c6 100644 --- a/packages/geotiff/README.md +++ b/packages/geotiff/README.md @@ -34,7 +34,7 @@ Until you try to load an image compressed with, say, [LERC], you don't pay for t [LERC]: https://github.com/Esri/lerc -### Transparent caching and chunking +### Full user control over caching and chunking There are a lot of great utilities in [`chunkd`](https://github.com/blacha/chunkd) that work out of the box here. diff --git a/packages/geotiff/package.json b/packages/geotiff/package.json index ca195c0..6533963 100644 --- a/packages/geotiff/package.json +++ b/packages/geotiff/package.json @@ -46,6 +46,10 @@ "extends": "../../package.json" }, "dependencies": { + "@chunkd/middleware": "^11.1.0", + "@chunkd/source": "^11.1.0", + "@chunkd/source-http": "^11.1.1", + "@chunkd/source-memory": "^11.0.2", "@cogeotiff/core": "^9.1.2", "@developmentseed/affine": "workspace:^", "@developmentseed/morecantile": "workspace:^", diff --git a/packages/geotiff/src/geotiff.ts b/packages/geotiff/src/geotiff.ts index 9a91051..f696441 100644 --- a/packages/geotiff/src/geotiff.ts +++ b/packages/geotiff/src/geotiff.ts @@ -1,11 +1,16 @@ +import { SourceCache } from "@chunkd/middleware"; +import { SourceChunk } from "@chunkd/middleware/build/src/middleware/chunk.js"; +import { SourceView } from "@chunkd/source"; +import { SourceHttp } from "@chunkd/source-http"; +import { SourceMemory } from "@chunkd/source-memory"; import type { Source, TiffImage } from "@cogeotiff/core"; import { Photometric, SubFileType, Tiff, TiffTag } from "@cogeotiff/core"; import type { Affine } from "@developmentseed/affine"; import type { ProjJson } from "./crs.js"; import { crsFromGeoKeys } from "./crs.js"; import { fetchTile } from "./fetch.js"; -import type { GeoKeyDirectory } from "./ifd.js"; -import { extractGeoKeyDirectory } from "./ifd.js"; +import type { CachedTags, GeoKeyDirectory } from "./ifd.js"; +import { extractGeoKeyDirectory, prefetchTags } from "./ifd.js"; import { Overview } from "./overview.js"; import type { Tile } from "./tile.js"; import { index, xy } from "./transform.js"; @@ -30,6 +35,9 @@ export class GeoTIFF { /** A cached CRS value. */ private _crs?: number | ProjJson; + /** Cached TIFF tags that are pre-fetched when opening the GeoTIFF. */ + readonly cachedTags: CachedTags; + /** The underlying Tiff instance. */ readonly tiff: Tiff; @@ -48,12 +56,14 @@ export class GeoTIFF { maskImage: TiffImage | null, gkd: GeoKeyDirectory, overviews: Overview[], + cachedTags: CachedTags, ) { this.tiff = tiff; this.image = image; this.maskImage = maskImage; this.gkd = gkd; this.overviews = overviews; + this.cachedTags = cachedTags; } /** @@ -115,9 +125,18 @@ export class GeoTIFF { return sb.width * sb.height - sa.width * sa.height; }); + const cachedTags = await prefetchTags(primaryImage); + // Two-phase construction: create the GeoTIFF first (with empty overviews), // then build Overviews that reference back to it. - const geotiff = new GeoTIFF(tiff, primaryImage, primaryMask, gkd, []); + const geotiff = new GeoTIFF( + tiff, + primaryImage, + primaryMask, + gkd, + [], + cachedTags, + ); const overviews: Overview[] = dataEntries.map(([key, dataImage]) => { const maskImage = maskIFDs.get(key) ?? null; @@ -130,6 +149,29 @@ export class GeoTIFF { return geotiff; } + static async fromArrayBuffer(input: ArrayBuffer): Promise { + const source = new SourceMemory("memory://input.tif", input); + return await GeoTIFF.create(source); + } + + static async fromUrl( + url: string | URL, + { + chunkSize = 32 * 1024, + cacheSize = 1024 * 1024 * 1024, + }: { chunkSize?: number; cacheSize?: number } = {}, + ): Promise { + // read files in chunks + const chunk = new SourceChunk({ size: chunkSize }); + // 1MB cache for recently accessed chunks + const cache = new SourceCache({ size: cacheSize }); + + const source = new SourceHttp(url); + const view = new SourceView(source, [chunk, cache]); + + return await GeoTIFF.create(view); + } + // ── Properties from the primary image ───────────────────────────────── /** diff --git a/packages/geotiff/src/ifd.ts b/packages/geotiff/src/ifd.ts index 13328f0..48127d6 100644 --- a/packages/geotiff/src/ifd.ts +++ b/packages/geotiff/src/ifd.ts @@ -1,21 +1,60 @@ -import type { - TiffImage, - TiffTag, - TiffTagGeoType, - TiffTagType, -} from "@cogeotiff/core"; -import { TiffTagGeo } from "@cogeotiff/core"; - -/** Subset of TIFF tags that are pre-fetched in {@link TiffImage.init}. */ -export interface PreFetchedTags { +import type { TiffImage, TiffTagGeoType, TiffTagType } from "@cogeotiff/core"; +import { SampleFormat, TiffTag, TiffTagGeo } from "@cogeotiff/core"; + +/** Subset of TIFF tags that we pre-fetch for easier visualization. */ +export interface CachedTags { + bitsPerSample: Uint16Array; + colorMap?: Uint16Array; // TiffTagType[TiffTag.ColorMap]; compression: TiffTagType[TiffTag.Compression]; - imageHeight: TiffTagType[TiffTag.ImageHeight]; - imageWidth: TiffTagType[TiffTag.ImageWidth]; - modelPixelScale?: TiffTagType[TiffTag.ModelPixelScale]; - modelTiePoint?: TiffTagType[TiffTag.ModelTiePoint]; - modelTransformation?: TiffTagType[TiffTag.ModelTransformation]; - tileHeight?: TiffTagType[TiffTag.TileHeight]; - tileWidth?: TiffTagType[TiffTag.TileWidth]; + nodata: number | null; + photometric: TiffTagType[TiffTag.Photometric]; + sampleFormat: TiffTagType[TiffTag.SampleFormat]; + samplesPerPixel: number; // TiffTagType[TiffTag.SamplesPerPixel]; +} + +/** Pre-fetch TIFF tags for easier visualization. */ +export async function prefetchTags(image: TiffImage): Promise { + // Compression is pre-fetched in init + const compression = image.value(TiffTag.Compression); + if (compression === null) { + throw new Error("Compression tag should always exist."); + } + + const nodata = image.noData; + + const [bitsPerSample, colorMap, photometric, sampleFormat, samplesPerPixel] = + await Promise.all([ + image.fetch(TiffTag.BitsPerSample), + image.fetch(TiffTag.ColorMap), + image.fetch(TiffTag.Photometric), + image.fetch(TiffTag.SampleFormat), + image.fetch(TiffTag.SamplesPerPixel), + ]); + + if (bitsPerSample === null) { + throw new Error("BitsPerSample tag should always exist."); + } + + if (samplesPerPixel === null) { + throw new Error("SamplesPerPixel tag should always exist."); + } + + if (photometric === null) { + throw new Error("Photometric tag should always exist."); + } + + return { + bitsPerSample: new Uint16Array(bitsPerSample), + colorMap: colorMap ? new Uint16Array(colorMap as number[]) : undefined, + compression, + nodata, + photometric, + // Uint is the default sample format according to the spec + // https://web.archive.org/web/20240329145340/https://www.awaresystems.be/imaging/tiff/tifftags/sampleformat.html + sampleFormat: sampleFormat ?? [SampleFormat.Uint], + // Waiting for release with https://github.com/blacha/cogeotiff/pull/1394 + samplesPerPixel: samplesPerPixel as number, + }; } /** diff --git a/packages/geotiff/src/tile-matrix-set.ts b/packages/geotiff/src/tile-matrix-set.ts index 9e65fd0..8c140ee 100644 --- a/packages/geotiff/src/tile-matrix-set.ts +++ b/packages/geotiff/src/tile-matrix-set.ts @@ -90,6 +90,18 @@ export function generateTileMatrixSet( const bbox = geotiff.bbox; const tr = geotiff.transform; + // Full-resolution level is appended last. + if (!geotiff.isTiled) { + throw new Error("GeoTIFF must be tiled to generate a TMS."); + } + + if (tr[1] !== 0 || tr[3] !== 0) { + // TileMatrixSet assumes orthogonal axes + throw new Error( + "COG TileMatrixSet with rotation/skewed geotransform is not supported", + ); + } + // Perhaps we should allow metersPerUnit to take any string const crsUnit = crs.units as | "m" @@ -131,11 +143,6 @@ export function generateTileMatrixSet( ); } - // Full-resolution level is appended last. - if (!geotiff.isTiled) { - throw new Error("GeoTIFF must be tiled to generate a TMS."); - } - tileMatrices.push( buildTileMatrix( String(geotiff.overviews.length), diff --git a/packages/geotiff/tests/geotiff.test.ts b/packages/geotiff/tests/geotiff.test.ts deleted file mode 100644 index 94dee3b..0000000 --- a/packages/geotiff/tests/geotiff.test.ts +++ /dev/null @@ -1,205 +0,0 @@ -import { Photometric, SubFileType } from "@cogeotiff/core"; -import { describe, expect, it } from "vitest"; -import { GeoTIFF, isMaskIfd } from "../src/geotiff.js"; -import { mockImage, mockTiff } from "./helpers.js"; - -describe("isMaskIfd", () => { - it("returns true for a mask IFD", () => { - const image = mockImage({ - width: 256, - height: 256, - subFileType: SubFileType.Mask, - photometric: Photometric.Mask, - }); - expect(isMaskIfd(image)).toBe(true); - }); - - it("returns false when SubFileType has no mask bit", () => { - const image = mockImage({ - width: 256, - height: 256, - subFileType: SubFileType.ReducedImage, - photometric: Photometric.Mask, - }); - expect(isMaskIfd(image)).toBe(false); - }); - - it("returns false when Photometric is not Mask", () => { - const image = mockImage({ - width: 256, - height: 256, - subFileType: SubFileType.Mask, - photometric: Photometric.MinIsBlack, - }); - expect(isMaskIfd(image)).toBe(false); - }); - - it("returns true when SubFileType combines ReducedImage + Mask bits", () => { - const image = mockImage({ - width: 256, - height: 256, - subFileType: SubFileType.ReducedImage | SubFileType.Mask, - photometric: Photometric.Mask, - }); - expect(isMaskIfd(image)).toBe(true); - }); - - it("returns false when SubFileType is absent (defaults to 0)", () => { - const image = mockImage({ - width: 256, - height: 256, - photometric: Photometric.Mask, - }); - expect(isMaskIfd(image)).toBe(false); - }); -}); - -describe("GeoTIFF", () => { - it("throws for empty TIFF", async () => { - const tiff = mockTiff([]); - await expect(GeoTIFF.fromTiff(tiff)).rejects.toThrow(/does not contain/); - }); - - it("creates a GeoTIFF from a single-image TIFF", async () => { - const primary = mockImage({ - width: 1000, - height: 1000, - origin: [0, 0, 0], - resolution: [1, -1, 0], - samplesPerPixel: 3, - }); - const tiff = mockTiff([primary]); - const geo = await GeoTIFF.fromTiff(tiff); - - expect(geo.width).toBe(1000); - expect(geo.height).toBe(1000); - expect(geo.count).toBe(3); - expect(geo.overviews).toHaveLength(0); - expect(geo.transform).toEqual([1, 0, 0, 0, -1, 0]); - }); - - it("classifies reduced-resolution IFDs as overviews", async () => { - const primary = mockImage({ - width: 1000, - height: 1000, - origin: [0, 0, 0], - resolution: [1, -1, 0], - }); - const ov1 = mockImage({ width: 500, height: 500 }); - const ov2 = mockImage({ width: 250, height: 250 }); - - const tiff = mockTiff([primary, ov1, ov2]); - const geo = await GeoTIFF.fromTiff(tiff); - - expect(geo.overviews).toHaveLength(2); - }); - - it("sorts overviews finest-to-coarsest", async () => { - const primary = mockImage({ - width: 1000, - height: 1000, - origin: [0, 0, 0], - resolution: [1, -1, 0], - }); - // Insert in reverse order - const small = mockImage({ width: 125, height: 125 }); - const medium = mockImage({ width: 250, height: 250 }); - const large = mockImage({ width: 500, height: 500 }); - - const tiff = mockTiff([primary, small, medium, large]); - const geo = await GeoTIFF.fromTiff(tiff); - - expect(geo.overviews).toHaveLength(3); - expect(geo.overviews[0]!.width).toBe(500); - expect(geo.overviews[1]!.width).toBe(250); - expect(geo.overviews[2]!.width).toBe(125); - }); - - it("separates mask IFDs from data IFDs", async () => { - const primary = mockImage({ - width: 1000, - height: 1000, - origin: [0, 0, 0], - resolution: [1, -1, 0], - }); - const ov = mockImage({ width: 500, height: 500 }); - const primaryMask = mockImage({ - width: 1000, - height: 1000, - subFileType: SubFileType.Mask, - photometric: Photometric.Mask, - }); - const ovMask = mockImage({ - width: 500, - height: 500, - subFileType: SubFileType.ReducedImage | SubFileType.Mask, - photometric: Photometric.Mask, - }); - - const tiff = mockTiff([primary, ov, primaryMask, ovMask]); - const geo = await GeoTIFF.fromTiff(tiff); - - // Only one data overview (the mask IFDs are paired, not listed as overviews) - expect(geo.overviews).toHaveLength(1); - expect(geo.overviews[0]!.maskImage).not.toBeNull(); - }); - - it("scales overview transforms correctly", async () => { - const primary = mockImage({ - width: 1000, - height: 1000, - origin: [100, 200, 0], - resolution: [0.01, -0.01, 0], - }); - const ov = mockImage({ width: 500, height: 500 }); - - const tiff = mockTiff([primary, ov]); - const geo = await GeoTIFF.fromTiff(tiff); - - const ovTransform = geo.overviews[0]!.transform; - // scale = 1000 / 500 = 2 - expect(ovTransform[0]).toBeCloseTo(0.02); // a * 2 - expect(ovTransform[4]).toBeCloseTo(-0.02); // e * 2 - // Origin unchanged - expect(ovTransform[2]).toBe(100); // c - expect(ovTransform[5]).toBe(200); // f - }); - - it("exposes nodata", async () => { - const primary = mockImage({ - width: 100, - height: 100, - origin: [0, 0, 0], - resolution: [1, -1, 0], - noData: -9999, - }); - const tiff = mockTiff([primary]); - const geo = await GeoTIFF.fromTiff(tiff); - expect(geo.nodata).toBe(-9999); - }); - - it("exposes bbox", async () => { - const primary = mockImage({ - width: 100, - height: 100, - origin: [0, 0, 0], - resolution: [1, -1, 0], - bbox: [-180, -90, 180, 90], - }); - const tiff = mockTiff([primary]); - const geo = await GeoTIFF.fromTiff(tiff); - expect(geo.bbox).toEqual([-180, -90, 180, 90]); - }); - - it("defaults count to 1 when SamplesPerPixel is absent", async () => { - const primary = mockImage({ - width: 100, - height: 100, - origin: [0, 0, 0], - resolution: [1, -1, 0], - }); - const tiff = mockTiff([primary]); - const geo = await GeoTIFF.fromTiff(tiff); - expect(geo.count).toBe(1); - }); -}); diff --git a/packages/geotiff/tests/helpers.ts b/packages/geotiff/tests/helpers.ts index ebf0b4a..1ffbfb4 100644 --- a/packages/geotiff/tests/helpers.ts +++ b/packages/geotiff/tests/helpers.ts @@ -1,7 +1,5 @@ import { resolve } from "node:path"; import { SourceFile } from "@chunkd/source-file"; -import type { Tiff, TiffImage } from "@cogeotiff/core"; -import { SampleFormat, TiffTag } from "@cogeotiff/core"; import { GeoTIFF } from "../src/geotiff.js"; // ── Fixture helpers ───────────────────────────────────────────────────── @@ -32,83 +30,3 @@ export async function loadGeoTIFF( const source = new SourceFile(path); return GeoTIFF.create(source); } - -// ── Mock helpers ──────────────────────────────────────────────────────── - -/** Create a mock TiffImage with configurable properties. */ -export function mockImage(opts: { - width: number; - height: number; - tileWidth?: number; - tileHeight?: number; - tiled?: boolean; - origin?: [number, number, number]; - resolution?: [number, number, number]; - subFileType?: number; - photometric?: number; - samplesPerPixel?: number; - noData?: number | null; - epsg?: number | null; - bbox?: [number, number, number, number]; - modelTransformation?: number[] | null; -}): TiffImage { - const tiled = opts.tiled ?? true; - const tags = new Map(); - - if (opts.subFileType != null) { - tags.set(TiffTag.SubFileType, opts.subFileType); - } - if (opts.photometric != null) { - tags.set(TiffTag.Photometric, opts.photometric); - } - if (opts.samplesPerPixel != null) { - tags.set(TiffTag.SamplesPerPixel, opts.samplesPerPixel); - } - if (opts.modelTransformation != null) { - tags.set(TiffTag.ModelTransformation, opts.modelTransformation); - } - - // Default SampleFormat and BitsPerSample so fetchTile works - if (!tags.has(TiffTag.SampleFormat)) { - tags.set(TiffTag.SampleFormat, [SampleFormat.Uint]); - } - if (!tags.has(TiffTag.BitsPerSample)) { - tags.set(TiffTag.BitsPerSample, [8]); - } - - return { - size: { width: opts.width, height: opts.height }, - tileSize: { - width: opts.tileWidth ?? 256, - height: opts.tileHeight ?? 256, - }, - isTiled: () => tiled, - origin: opts.origin ?? [0, 0, 0], - resolution: opts.resolution ?? [1, -1, 0], - noData: opts.noData ?? null, - epsg: opts.epsg ?? null, - bbox: opts.bbox ?? [0, 0, 100, 100], - init: async () => {}, - has: (tag: number) => tags.has(tag), - value: (tag: number) => { - if (tags.has(tag)) return tags.get(tag); - return null; - }, - valueGeo: () => null, - fetch: async (tag: number) => tags.get(tag) ?? null, - getTile: async (_x: number, _y: number) => ({ - bytes: new ArrayBuffer( - (opts.tileWidth ?? 256) * - (opts.tileHeight ?? 256) * - (opts.samplesPerPixel ?? 1), - ), - mimeType: "application/octet-stream", - compression: 1, // None - }), - } as unknown as TiffImage; -} - -/** Create a mock Tiff with the given images. */ -export function mockTiff(images: TiffImage[]): Tiff { - return { images } as unknown as Tiff; -} diff --git a/packages/geotiff/tests/overview.test.ts b/packages/geotiff/tests/overview.test.ts deleted file mode 100644 index 170e50f..0000000 --- a/packages/geotiff/tests/overview.test.ts +++ /dev/null @@ -1,26 +0,0 @@ -import { describe, expect, it } from "vitest"; -import { GeoTIFF } from "../src/geotiff.js"; -import { mockImage, mockTiff } from "./helpers.js"; - -describe("Overview", () => { - it("computes scaled transform from parent", async () => { - const primary = mockImage({ - width: 1000, - height: 1000, - origin: [100, 200, 0], - resolution: [0.01, -0.01, 0], - }); - const ov = mockImage({ width: 500, height: 500 }); - - const tiff = mockTiff([primary, ov]); - const geo = await GeoTIFF.fromTiff(tiff); - - const ovTransform = geo.overviews[0]!.transform; - // scale = 1000 / 500 = 2 - expect(ovTransform[0]).toBeCloseTo(0.02); // a * 2 - expect(ovTransform[4]).toBeCloseTo(-0.02); // e * 2 - // Origin unchanged - expect(ovTransform[2]).toBe(100); // c - expect(ovTransform[5]).toBe(200); // f - }); -}); diff --git a/packages/geotiff/tests/tile-matrix-set.test.ts b/packages/geotiff/tests/tile-matrix-set.test.ts index 5cdfdb2..be5fb63 100644 --- a/packages/geotiff/tests/tile-matrix-set.test.ts +++ b/packages/geotiff/tests/tile-matrix-set.test.ts @@ -1,5 +1,6 @@ import { describe, expect, it } from "vitest"; import wktParser from "wkt-parser"; +import { GeoTIFF } from "../src/geotiff.js"; import { generateTileMatrixSet } from "../src/tile-matrix-set.js"; import { loadGeoTIFF } from "./helpers.js"; @@ -126,3 +127,85 @@ describe("test TMS", () => { ]); }); }); + +describe("create TileMatrixSet from COG", () => { + it("creates TMS", async () => { + const url = + "https://ds-wheels.s3.us-east-1.amazonaws.com/m_4007307_sw_18_060_20220803.tif"; + + const geotiff = await GeoTIFF.fromUrl(url); + + const tms = generateTileMatrixSet(geotiff, { units: "m" }); + + const expectedTileMatrices = [ + { + id: "0", + scaleDenominator: 68684.95742667928, + cellSize: 19.231788079470196, + pointOfOrigin: [647118, 4533600], + tileWidth: 512, + tileHeight: 512, + matrixWidth: 1, + matrixHeight: 1, + cornerOfOrigin: "topLeft", + }, + { + id: "1", + scaleDenominator: 34285.71428571429, + cellSize: 9.6, + pointOfOrigin: [647118, 4533600], + tileWidth: 512, + tileHeight: 512, + matrixWidth: 2, + matrixHeight: 2, + cornerOfOrigin: "topLeft", + }, + { + id: "2", + scaleDenominator: 17142.857142857145, + cellSize: 4.8, + pointOfOrigin: [647118, 4533600], + tileWidth: 512, + tileHeight: 512, + matrixWidth: 3, + matrixHeight: 4, + cornerOfOrigin: "topLeft", + }, + { + id: "3", + scaleDenominator: 8571.428571428572, + cellSize: 2.4, + pointOfOrigin: [647118, 4533600], + tileWidth: 512, + tileHeight: 512, + matrixWidth: 5, + matrixHeight: 7, + cornerOfOrigin: "topLeft", + }, + { + id: "4", + scaleDenominator: 4285.714285714286, + cellSize: 1.2, + pointOfOrigin: [647118, 4533600], + tileWidth: 512, + tileHeight: 512, + matrixWidth: 10, + matrixHeight: 13, + cornerOfOrigin: "topLeft", + }, + { + id: "5", + scaleDenominator: 2142.857142857143, + cellSize: 0.6, + pointOfOrigin: [647118, 4533600], + tileWidth: 512, + tileHeight: 512, + matrixWidth: 19, + matrixHeight: 25, + cornerOfOrigin: "topLeft", + }, + ]; + + expect(tms.tileMatrices).toStrictEqual(expectedTileMatrices); + }); +}); diff --git a/packages/geotiff/tsconfig.build.json b/packages/geotiff/tsconfig.build.json index 337c3a4..0fb6fc5 100644 --- a/packages/geotiff/tsconfig.build.json +++ b/packages/geotiff/tsconfig.build.json @@ -7,5 +7,8 @@ }, "include": ["src/**/*"], "exclude": ["node_modules", "dist", "tests"], - "references": [{ "path": "../affine/tsconfig.build.json" }] + "references": [ + { "path": "../affine/tsconfig.build.json" }, + { "path": "../morecantile/tsconfig.build.json" } + ] } diff --git a/packages/morecantile/package.json b/packages/morecantile/package.json index a82e6e3..3004203 100644 --- a/packages/morecantile/package.json +++ b/packages/morecantile/package.json @@ -45,7 +45,9 @@ "typescript": "^5.9.3", "vitest": "^4.0.18" }, - "dependencies": {}, + "dependencies": { + "@developmentseed/affine": "workspace:^" + }, "peerDependencies": {}, "volta": { "extends": "../../package.json" diff --git a/packages/morecantile/src/index.ts b/packages/morecantile/src/index.ts index 155e30d..c48db5d 100644 --- a/packages/morecantile/src/index.ts +++ b/packages/morecantile/src/index.ts @@ -1,4 +1,4 @@ -export type { Affine } from "./transform.js"; +export { xy_bounds } from "./tile.js"; export { matrixTransform, tileTransform } from "./transform.js"; export type { BoundingBox, diff --git a/packages/morecantile/src/tile.ts b/packages/morecantile/src/tile.ts new file mode 100644 index 0000000..f9d16af --- /dev/null +++ b/packages/morecantile/src/tile.ts @@ -0,0 +1,64 @@ +import * as affine from "@developmentseed/affine"; +import { tileTransform } from "./transform"; +import type { BoundingBox, TileMatrix, TileMatrixSet } from "./types"; +import { narrowTileMatrixSet } from "./utils"; + +/** + * Return the bounding box of the tile in the TMS's native coordinate reference + * system. + */ +export function xy_bounds( + matrix: TileMatrix, + tile: { x: number; y: number }, +): BoundingBox; +/** + * Return the bounding box of the tile in the TMS's native coordinate reference + * system. + */ +export function xy_bounds( + matrixSet: TileMatrixSet, + tile: { x: number; y: number; z: number }, +): BoundingBox; +export function xy_bounds( + input: TileMatrix | TileMatrixSet, + tile: { x: number; y: number; z?: number }, +): BoundingBox { + const tileMatrix = getTileMatrix(input, tile); + const { tileHeight, tileWidth } = tileMatrix; + const { x, y } = tile; + const tileAffine = tileTransform(tileMatrix, { col: x, row: y }); + + // Apply affine to local tile pixel corners (0,0) is the origin corner, + // (tileWidth, tileHeight) is the opposite corner. + const [x0, y0] = affine.apply(tileAffine, 0, 0); + const [x1, y1] = affine.apply(tileAffine, tileWidth, tileHeight); + + if (tileMatrix.cornerOfOrigin === "bottomLeft") { + // (x0, y0) is bottom-left, (x1, y1) is top-right + return { lowerLeft: [x0, y0], upperRight: [x1, y1] }; + } + + // topLeft (default): (x0, y0) is top-left, (x1, y1) is bottom-right + return { lowerLeft: [x0, y1], upperRight: [x1, y0] }; +} + +function getTileMatrix( + input: TileMatrix | TileMatrixSet, + tile: { x: number; y: number; z?: number }, +): TileMatrix { + if (narrowTileMatrixSet(input)) { + if (tile.z === undefined) { + throw new Error("Tile z level is required when input is a TileMatrixSet"); + } + const tileMatrix = input.tileMatrices[tile.z]; + if (!tileMatrix) { + throw new Error( + `Tile z level ${tile.z} is out of bounds for TileMatrixSet with ${input.tileMatrices.length} levels.`, + ); + } + + return tileMatrix; + } else { + return input; + } +} diff --git a/packages/morecantile/src/transform.ts b/packages/morecantile/src/transform.ts index 5ce94f5..1725dcf 100644 --- a/packages/morecantile/src/transform.ts +++ b/packages/morecantile/src/transform.ts @@ -1,28 +1,6 @@ +import type { Affine } from "@developmentseed/affine"; import type { TileMatrix } from "./types/index.js"; -/** - * A 2-D affine transform as a six-element tuple in row-major order: - * - * x_crs = a * col_px + b * row_px + c - * y_crs = d * col_px + e * row_px + f - * - * For the north-up, axis-aligned grids produced by OGC TileMatrices, - * b and d are always zero, but the type does not enforce that. - * - * Note: the meaning of the two output axes depends on the - * `orderedAxes` field of the parent TileMatrixSet, which these - * helpers do not consult. Axis 0 of pointOfOrigin is axis 0 of the - * transform output, whatever that axis happens to be. - */ -export type Affine = [ - a: number, - b: number, - c: number, - d: number, - e: number, - f: number, -]; - /** * Construct a single affine transform that maps pixel coordinates * within *any* tile of the matrix to CRS coordinates. diff --git a/packages/morecantile/src/utils.ts b/packages/morecantile/src/utils.ts index 4a00f1e..5e478ed 100644 --- a/packages/morecantile/src/utils.ts +++ b/packages/morecantile/src/utils.ts @@ -12,6 +12,9 @@ * @param semiMajorAxis - The semi-major axis of the ellipsoid, required if unit is 'degree'. * @returns The meters per unit conversion factor. */ + +import type { TileMatrix, TileMatrixSet } from "./types"; + // https://github.com/developmentseed/morecantile/blob/7c95a11c491303700d6e33e9c1607f2719584dec/morecantile/utils.py#L67-L90 export function metersPerUnit( unit: @@ -22,7 +25,7 @@ export function metersPerUnit( | "foot" | "us survey foot" | "degree", - { semiMajorAxis }: { semiMajorAxis?: number }, + { semiMajorAxis }: { semiMajorAxis?: number } = {}, ): number { unit = unit.toLowerCase() as typeof unit; switch (unit) { @@ -52,3 +55,9 @@ export function metersPerUnit( `Unsupported CRS units: ${unit} when computing metersPerUnit.`, ); } + +export function narrowTileMatrixSet( + obj: TileMatrix | TileMatrixSet, +): obj is TileMatrixSet { + return "tileMatrices" in obj; +} diff --git a/packages/morecantile/tests/utils.test.ts b/packages/morecantile/tests/utils.test.ts new file mode 100644 index 0000000..06d392a --- /dev/null +++ b/packages/morecantile/tests/utils.test.ts @@ -0,0 +1,13 @@ +import { describe, expect, it } from "vitest"; +import { metersPerUnit } from "../src"; + +describe("metersPerUnit", () => { + it("handles lowercase us survey foot", () => { + expect(metersPerUnit("us survey foot")).toBe(1200 / 3937); + }); + + it("handles mixed case US Survey Foot", () => { + // @ts-expect-error testing case insensitivity with non-standard casing + expect(metersPerUnit("US Survey Foot")).toBe(1200 / 3937); + }); +}); diff --git a/packages/raster-reproject/package.json b/packages/raster-reproject/package.json index a7901bd..120c47e 100644 --- a/packages/raster-reproject/package.json +++ b/packages/raster-reproject/package.json @@ -9,10 +9,6 @@ ".": { "types": "./dist/index.d.ts", "default": "./dist/index.js" - }, - "./affine": { - "types": "./dist/affine/index.d.ts", - "default": "./dist/affine/index.js" } }, "sideEffects": false, @@ -44,6 +40,7 @@ "url": "git+https://github.com/developmentseed/deck.gl-raster.git" }, "devDependencies": { + "@developmentseed/affine": "workspace:^", "@types/node": "^25.1.0", "jsdom": "^27.4.0", "proj4": "^2.20.2", diff --git a/packages/raster-reproject/src/affine/affine.ts b/packages/raster-reproject/src/affine/affine.ts deleted file mode 100644 index d455021..0000000 --- a/packages/raster-reproject/src/affine/affine.ts +++ /dev/null @@ -1,51 +0,0 @@ -export type GeoTransform = [number, number, number, number, number, number]; - -/** - * Find the inverse of this GeoTransform. - * - * Ported from rasterio/affine: - * https://github.com/rasterio/affine/blob/a7a916fc7012f8afeb6489246ada61a76ccb8bc7/src/affine.py#L671-L692 - * under the BSD-3-Clause License. - * - * @param {GeoTransform} gt Geotransform. - * - * @return {GeoTransform} Inverse of the geotransform. - */ -export function invertGeoTransform(gt: GeoTransform): GeoTransform { - if (isDegenerate(gt)) { - throw new Error("Cannot invert degenerate transform"); - } - - const idet = 1.0 / determinant(gt); - const [sa, sb, sc, sd, se, sf] = gt; - const ra = se * idet; - const rb = -sb * idet; - const rd = -sd * idet; - const re = sa * idet; - // biome-ignore format: array - return [ - ra, rb, -sc * ra - sf * rb, - rd, re, -sc * rd - sf * re, - ]; -} - -function isDegenerate(gt: GeoTransform): boolean { - return determinant(gt) === 0; -} - -function determinant(gt: GeoTransform): number { - const [a, b, _c, d, e, _f] = gt; - return a * e - b * d; -} - -/** - * Apply a GeoTransform to a coordinate. - */ -export function applyAffine( - x: number, - y: number, - gt: [number, number, number, number, number, number], -): [number, number] { - const [a, b, c, d, e, f] = gt; - return [a * x + b * y + c, d * x + e * y + f]; -} diff --git a/packages/raster-reproject/src/affine/index.ts b/packages/raster-reproject/src/affine/index.ts deleted file mode 100644 index 7a082af..0000000 --- a/packages/raster-reproject/src/affine/index.ts +++ /dev/null @@ -1,2 +0,0 @@ -export type { GeoTransform } from "./affine.js"; -export { applyAffine, invertGeoTransform } from "./affine.js"; diff --git a/packages/raster-reproject/tests/example.test.ts b/packages/raster-reproject/tests/example.test.ts index ee28e05..db689ae 100644 --- a/packages/raster-reproject/tests/example.test.ts +++ b/packages/raster-reproject/tests/example.test.ts @@ -1,11 +1,8 @@ import { readFileSync, writeFileSync } from "node:fs"; import { dirname, join, resolve } from "node:path"; import { fileURLToPath } from "node:url"; +import * as affine from "@developmentseed/affine"; import { RasterReprojector } from "@developmentseed/raster-reproject"; -import { - applyAffine, - invertGeoTransform, -} from "@developmentseed/raster-reproject/affine"; import proj4 from "proj4"; import type { PROJJSONDefinition } from "proj4/dist/lib/core"; import { describe, it } from "vitest"; @@ -25,17 +22,18 @@ type FixtureJSON = { }; // Note: this is copied from deck.gl-geotiff -function fromGeoTransform( +function fromAffine( geotransform: [number, number, number, number, number, number], ): { forwardTransform: (x: number, y: number) => [number, number]; inverseTransform: (x: number, y: number) => [number, number]; } { - const inverseGeotransform = invertGeoTransform(geotransform); + const inverseGeotransform = affine.invert(geotransform); return { - forwardTransform: (x: number, y: number) => applyAffine(x, y, geotransform), + forwardTransform: (x: number, y: number) => + affine.apply(geotransform, x, y), inverseTransform: (x: number, y: number) => - applyAffine(x, y, inverseGeotransform), + affine.apply(inverseGeotransform, x, y), }; } @@ -66,8 +64,7 @@ function parseFixture(fixturePath: string): RasterReprojector { affineGeotransform = geotransform; } - const { inverseTransform, forwardTransform } = - fromGeoTransform(affineGeotransform); + const { inverseTransform, forwardTransform } = fromAffine(affineGeotransform); const converter = proj4(projjson || wkt2, "EPSG:4326"); const reprojectionFns = { diff --git a/pnpm-lock.yaml b/pnpm-lock.yaml index f4755f1..59d1280 100644 --- a/pnpm-lock.yaml +++ b/pnpm-lock.yaml @@ -57,18 +57,15 @@ importers: '@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) - geotiff: - specifier: ^2.1.3 - version: 2.1.3 - geotiff-geokeys-to-proj4: - specifier: ^2024.4.13 - version: 2024.4.13 maplibre-gl: specifier: ^5.17.0 version: 5.17.0 @@ -243,6 +240,9 @@ importers: packages/deck.gl-geotiff: dependencies: + '@cogeotiff/core': + specifier: ^9.1.2 + version: 9.1.2 '@deck.gl/core': specifier: 9.2.7 version: 9.2.7 @@ -255,9 +255,18 @@ importers: '@deck.gl/mesh-layers': specifier: 9.2.7 version: 9.2.7(@deck.gl/core@9.2.7)(@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)) + '@developmentseed/affine': + specifier: workspace:^ + version: link:../affine '@developmentseed/deck.gl-raster': specifier: workspace:^ version: link:../deck.gl-raster + '@developmentseed/geotiff': + specifier: workspace:^ + version: link:../geotiff + '@developmentseed/morecantile': + specifier: workspace:^ + version: link:../morecantile '@developmentseed/raster-reproject': specifier: workspace:^ version: link:../raster-reproject @@ -267,12 +276,12 @@ importers: flatbush: specifier: ^4.5.0 version: 4.5.0 - geotiff: - specifier: 2.1.3 - version: 2.1.3 proj4: specifier: ^2.20.2 version: 2.20.2 + wkt-parser: + specifier: ^1.5.2 + version: 1.5.2 devDependencies: '@types/node': specifier: ^25.1.0 @@ -301,6 +310,12 @@ importers: '@deck.gl/mesh-layers': specifier: 9.2.7 version: 9.2.7(@deck.gl/core@9.2.7)(@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)) + '@developmentseed/affine': + specifier: workspace:^ + version: link:../affine + '@developmentseed/morecantile': + specifier: workspace:^ + version: link:../morecantile '@developmentseed/raster-reproject': specifier: workspace:^ version: link:../raster-reproject @@ -363,6 +378,18 @@ importers: packages/geotiff: dependencies: + '@chunkd/middleware': + specifier: ^11.1.0 + version: 11.1.0 + '@chunkd/source': + specifier: ^11.1.0 + version: 11.1.0 + '@chunkd/source-http': + specifier: ^11.1.1 + version: 11.1.1 + '@chunkd/source-memory': + specifier: ^11.0.2 + version: 11.0.2 '@cogeotiff/core': specifier: ^9.1.2 version: 9.1.2 @@ -396,6 +423,10 @@ importers: version: 1.5.2 packages/morecantile: + dependencies: + '@developmentseed/affine': + specifier: workspace:^ + version: link:../affine devDependencies: '@types/node': specifier: ^25.1.0 @@ -418,6 +449,9 @@ importers: packages/raster-reproject: devDependencies: + '@developmentseed/affine': + specifier: workspace:^ + version: link:../affine '@types/node': specifier: ^25.1.0 version: 25.1.0 @@ -588,10 +622,22 @@ packages: cpu: [x64] os: [win32] + '@chunkd/middleware@11.1.0': + resolution: {integrity: sha512-xLEoDWDG3djf/plcsEsm0q6eAD3fFC/nJ6QV9234a+DE4iN8xuaqdHqtk9CRC877YfCpWZ1RDw3T62OUeqyWqQ==} + engines: {node: '>=16.0.0'} + '@chunkd/source-file@11.0.1': resolution: {integrity: sha512-9DoA1djozuDpUklg0CRj/GBEJzvZeqwr1EIOTcEr3igDAf+UZ3y3lOgUaIWuNeATbwQ7Pdml2S6SLLAdAe45wQ==} engines: {node: '>=16.0.0'} + '@chunkd/source-http@11.1.1': + resolution: {integrity: sha512-cJPPYQ8JB9y+zWvv6NKOLdOp+jXYejTuzLa9bCtwFunlKML7a4PuOZVAQp8xlke6+HOq4w6Dzednq7mhgzHKaA==} + engines: {node: '>=18.0.0'} + + '@chunkd/source-memory@11.0.2': + resolution: {integrity: sha512-K2yJTvLv+Z49YKGz43Gb2Q1/6ARtnUUh9EgVgqKZS/Ce9ND9IXkgmfyrYShY2/d9iJ4B7Nq1LZIJubDBhoYtww==} + engines: {node: '>=16.0.0'} + '@chunkd/source@11.1.0': resolution: {integrity: sha512-3BcrHK4XDm3t4vDx1OisgcOZ05+EarEqwaGVIL7IMHUmkXwN/Aprlnr7aVP3BYcltgcCcfMd1pXQicbSrP563g==} engines: {node: '>=16.0.0'} @@ -2762,10 +2808,22 @@ snapshots: '@biomejs/cli-win32-x64@2.3.13': optional: true + '@chunkd/middleware@11.1.0': + dependencies: + '@chunkd/source': 11.1.0 + '@chunkd/source-file@11.0.1': dependencies: '@chunkd/source': 11.1.0 + '@chunkd/source-http@11.1.1': + dependencies: + '@chunkd/source': 11.1.0 + + '@chunkd/source-memory@11.0.2': + dependencies: + '@chunkd/source': 11.1.0 + '@chunkd/source@11.1.0': {} '@cogeotiff/core@9.1.2': {} diff --git a/tsconfig.json b/tsconfig.json index 2bbec57..aa1013e 100644 --- a/tsconfig.json +++ b/tsconfig.json @@ -27,11 +27,11 @@ "files": [], "references": [ { "path": "packages/affine/tsconfig.build.json" }, - { "path": "packages/morecantile/tsconfig.build.json" }, { "path": "packages/deck.gl-geotiff/tsconfig.build.json" }, { "path": "packages/deck.gl-raster/tsconfig.build.json" }, { "path": "packages/deck.gl-zarr/tsconfig.build.json" }, { "path": "packages/geotiff/tsconfig.build.json" }, + { "path": "packages/morecantile/tsconfig.build.json" }, { "path": "packages/raster-reproject/tsconfig.build.json" } ] }