From 87fa033fb7669f1c32269d56bc6772df44cc0212 Mon Sep 17 00:00:00 2001 From: Wentao Kuang Date: Mon, 13 Jan 2025 15:12:24 +1300 Subject: [PATCH] create topo raster cog --- packages/cogify/src/cogify/cli/cli.cog.ts | 95 +++++++++++++----- packages/cogify/src/cogify/gdal.command.ts | 110 ++++++++++++--------- packages/cogify/src/cogify/stac.ts | 35 +------ packages/cogify/src/tile.cover.ts | 1 - 4 files changed, 135 insertions(+), 106 deletions(-) diff --git a/packages/cogify/src/cogify/cli/cli.cog.ts b/packages/cogify/src/cogify/cli/cli.cog.ts index 41ea845e5..ec0b83ab1 100644 --- a/packages/cogify/src/cogify/cli/cli.cog.ts +++ b/packages/cogify/src/cogify/cli/cli.cog.ts @@ -1,5 +1,5 @@ import { isEmptyTiff } from '@basemaps/config-loader'; -import { Projection, ProjectionLoader, TileMatrixSet, TileMatrixSets } from '@basemaps/geo'; +import { Projection, ProjectionLoader, TileId, TileMatrixSet, TileMatrixSets } from '@basemaps/geo'; import { fsa, LogType, stringToUrlFolder, Tiff } from '@basemaps/shared'; import { CliId, CliInfo } from '@basemaps/shared/build/cli/info.js'; import { Metrics } from '@linzjs/metrics'; @@ -15,7 +15,13 @@ import { CutlineOptimizer } from '../../cutline.js'; import { SourceDownloader } from '../../download.js'; import { HashTransform } from '../../hash.stream.js'; import { getLogger, logArguments } from '../../log.js'; -import { gdalBuildCog, gdalBuildVrt, gdalBuildVrtWarp, gdalCreate } from '../gdal.command.js'; +import { + gdalBuildCog, + gdalBuildTopoRasterCommands, + gdalBuildVrt, + gdalBuildVrtWarp, + gdalCreate, +} from '../gdal.command.js'; import { GdalRunner } from '../gdal.runner.js'; import { Url, UrlArrayJsonFile } from '../parsers.js'; import { CogifyCreationOptions, CogifyStacItem, getCutline, getSources } from '../stac.js'; @@ -65,6 +71,7 @@ export const BasemapsCogifyCreateCommand = command({ ...logArguments, path: restPositionals({ type: Url, displayName: 'path', description: 'Path to covering configuration' }), force: flag({ long: 'force', description: 'Overwrite existing tiff files' }), + topo: flag({ long: 'topo', description: 'Standardize topo raster map' }), concurrency: option({ type: number, long: 'concurrency', @@ -155,7 +162,7 @@ export const BasemapsCogifyCreateCommand = command({ const { item, url } = f; const cutlineLink = getCutline(item.links); const options = item.properties['linz_basemaps:options']; - const tileId = options.tileId; + const tileId = args.topo ? item.id : TileId.fromTile(options.tile); // Location to where the tiff should be stored const tiffPath = new URL(tileId + '.tiff', url); @@ -173,11 +180,24 @@ export const BasemapsCogifyCreateCommand = command({ // Create the tiff concurrently const outputTiffPath = await Q(async () => { metrics.start(tileId); // Only start the timer when the cog is actually being processed - - const cutline = await CutlineOptimizer.loadFromLink(cutlineLink, tileMatrix); const sourceLocations = await Promise.all(sourceFiles.map((f) => sources.get(f, logger))); + const cutline = await CutlineOptimizer.loadFromLink(cutlineLink, tileMatrix); + if (args.topo) { + const width = item.properties['source.width'] as number; + const height = item.properties['source.height'] as number; + return createTopoCog({ + tileId, + options, + tempFolder: tmpFolder, + sourceFiles: sourceLocations, + cutline, + size: { width, height }, + logger, + }); + } return createCog({ + tileId, options, tempFolder: tmpFolder, sourceFiles: sourceLocations, @@ -268,7 +288,9 @@ export const BasemapsCogifyCreateCommand = command({ { count: toCreate.length, created: filtered.length, - files: filtered.map((f) => f.item.properties['linz_basemaps:options'].tileId), + files: filtered.map((f) => { + return args.topo ? f.item.id : TileId.fromTile(f.item.properties['linz_basemaps:options'].tile); + }), }, 'Cog:Done', ); @@ -276,6 +298,8 @@ export const BasemapsCogifyCreateCommand = command({ }); export interface CogCreationContext { + /** TileId for the file name */ + tileId: string; /** COG Creation options */ options: CogifyCreationOptions; /** Location to store all the temporary files */ @@ -284,6 +308,8 @@ export interface CogCreationContext { sourceFiles: URL[]; /** Optional cutline to cut the imagery too */ cutline: CutlineOptimizer; + /** Optional Source imagery size for topo raster trim pixel */ + size?: { width: number; height: number }; /** Optional logger */ logger?: LogType; } @@ -292,7 +318,7 @@ export interface CogCreationContext { async function createCog(ctx: CogCreationContext): Promise { const options = ctx.options; await ProjectionLoader.load(options.sourceEpsg); - const tileId = options.tileId; + const tileId = ctx.tileId; const logger = ctx.logger?.child({ tileId }); @@ -303,13 +329,13 @@ async function createCog(ctx: CogCreationContext): Promise { logger?.debug({ tileId }, 'Cog:Create:VrtSource'); // Create the vrt of all the source files - const vrtSourceCommand = gdalBuildVrt(new URL(`${tileId}-source.vrt`, ctx.tempFolder), ctx.sourceFiles, options); + const vrtSourceCommand = gdalBuildVrt(new URL(`${tileId}-source.vrt`, ctx.tempFolder), ctx.sourceFiles); await new GdalRunner(vrtSourceCommand).run(logger); logger?.debug({ tileId }, 'Cog:Create:VrtWarp'); const cutlineProperties: { url: URL | null; blend: number } = { url: null, blend: ctx.cutline.blend }; - if (ctx.cutline.path && options.tile) { + if (ctx.cutline) { logger?.debug('Cog:Cutline'); const optimizedCutline = ctx.cutline.optimize(options.tile); if (optimizedCutline) { @@ -321,23 +347,19 @@ async function createCog(ctx: CogCreationContext): Promise { } } - let vrtOutput = vrtSourceCommand.output; - if (!options.noReprojecting) { - // warp the source VRT into the output parameters - const vrtWarpCommand = gdalBuildVrtWarp( - new URL(`${tileId}-${options.tileMatrix}-warp.vrt`, ctx.tempFolder), - vrtSourceCommand.output, - options.sourceEpsg, - cutlineProperties, - options, - ); - await new GdalRunner(vrtWarpCommand).run(logger); - vrtOutput = vrtWarpCommand.output; - } + // warp the source VRT into the output parameters + const vrtWarpCommand = gdalBuildVrtWarp( + new URL(`${tileId}-${options.tileMatrix}-warp.vrt`, ctx.tempFolder), + vrtSourceCommand.output, + options.sourceEpsg, + cutlineProperties, + options, + ); + await new GdalRunner(vrtWarpCommand).run(logger); if (options.background == null) { // Create the COG from the warped vrt without a forced background - const cogCreateCommand = gdalBuildCog(new URL(`${tileId}.tiff`, ctx.tempFolder), vrtOutput, options); + const cogCreateCommand = gdalBuildCog(new URL(`${tileId}.tiff`, ctx.tempFolder), vrtWarpCommand.output, options); await new GdalRunner(cogCreateCommand).run(logger); return cogCreateCommand.output; } @@ -349,7 +371,7 @@ async function createCog(ctx: CogCreationContext): Promise { // Create a vrt with the background tiff behind the source file vrt const vrtMergeCommand = gdalBuildVrt(new URL(`${tileId}-merged.vrt`, ctx.tempFolder), [ gdalCreateCommand.output, - vrtOutput, + vrtWarpCommand.output, ]); await new GdalRunner(vrtMergeCommand).run(logger); @@ -359,6 +381,31 @@ async function createCog(ctx: CogCreationContext): Promise { return cogCreateCommand.output; } +/** Create a cog from the creation options */ +async function createTopoCog(ctx: CogCreationContext): Promise { + const options = ctx.options; + await ProjectionLoader.load(options.sourceEpsg); + const tileId = ctx.tileId; + + const logger = ctx.logger?.child({ tileId }); + + logger?.debug({ tileId }, 'TopoCog:Create:VrtSource'); + // Create the vrt of all the source files + const vrtSourceCommand = gdalBuildVrt(new URL(`${tileId}-source.vrt`, ctx.tempFolder), ctx.sourceFiles, true); + await new GdalRunner(vrtSourceCommand).run(logger); + + // Create the COG from the vrt file + if (ctx.size == null) throw new Error('TopoCog: Source image size is required for pixel trim'); + const cogCreateCommand = gdalBuildTopoRasterCommands( + new URL(`${tileId}.tiff`, ctx.tempFolder), + vrtSourceCommand.output, + ctx.size?.width, + ctx.size?.height, + ); + await new GdalRunner(cogCreateCommand).run(logger); + return cogCreateCommand.output; +} + /** * Very basic checking for the output tiff to ensure it was uploaded ok * Just open it as a COG and ensure the metadata looks about right diff --git a/packages/cogify/src/cogify/gdal.command.ts b/packages/cogify/src/cogify/gdal.command.ts index 9831c495f..8e37398ce 100644 --- a/packages/cogify/src/cogify/gdal.command.ts +++ b/packages/cogify/src/cogify/gdal.command.ts @@ -7,44 +7,14 @@ import { GdalCommand } from './gdal.runner.js'; import { CogifyCreationOptions } from './stac.js'; const isPowerOfTwo = (x: number): boolean => (x & (x - 1)) === 0; +const DEFAULT_TRIM_PIXEL_RIGHT = 1.7; // 1.7 pixels to trim from the right side of the topo raster imagery -interface TargetOptions { - targetSrs?: string; - extent?: number[]; - targetResolution?: number; -} - -function getTargetOptions(opt: CogifyCreationOptions): TargetOptions { - const targetOpts: TargetOptions = {}; - - if (opt.tileMatrix) { - const tileMatrix = TileMatrixSets.find(opt.tileMatrix); - if (tileMatrix == null) throw new Error('Unable to find tileMatrix: ' + opt.tileMatrix); - if (opt.noReprojecting == null) targetOpts.targetSrs = tileMatrix.projection.toEpsgString(); - - if (opt.tile) { - const bounds = tileMatrix.tileToSourceBounds(opt.tile); - targetOpts.extent = [ - Math.min(bounds.x, bounds.right), - Math.min(bounds.y, bounds.bottom), - Math.max(bounds.x, bounds.right), - Math.max(bounds.y, bounds.bottom), - ]; - } - - if (opt.zoomLevel) { - targetOpts.targetResolution = tileMatrix.pixelScale(opt.zoomLevel); - } - } - return targetOpts; -} - -export function gdalBuildVrt(targetVrt: URL, source: URL[], opt?: CogifyCreationOptions): GdalCommand { +export function gdalBuildVrt(targetVrt: URL, source: URL[], addalpha?: boolean): GdalCommand { if (source.length === 0) throw new Error('No source files given for :' + targetVrt.href); return { output: targetVrt, command: 'gdalbuildvrt', - args: [opt && opt.addalpha ? ['-addalpha'] : undefined, urlToString(targetVrt), ...source.map(urlToString)] + args: [addalpha ? ['-addalpha'] : undefined, urlToString(targetVrt), ...source.map(urlToString)] .filter((f) => f != null) .flat() .map(String), @@ -60,7 +30,6 @@ export function gdalBuildVrtWarp( ): GdalCommand { const tileMatrix = TileMatrixSets.find(opt.tileMatrix); if (tileMatrix == null) throw new Error('Unable to find tileMatrix: ' + opt.tileMatrix); - if (opt.zoomLevel == null) throw new Error('Unable to find zoomLevel'); const targetResolution = tileMatrix.pixelScale(opt.zoomLevel); return { @@ -68,6 +37,8 @@ export function gdalBuildVrtWarp( command: 'gdalwarp', args: [ ['-of', 'vrt'], // Output as a VRT + // ['-co', 'compress=lzw'], + // ['-co', 'bigtiff=yes'], '-multi', // Mutithread IO ['-wo', 'NUM_THREADS=ALL_CPUS'], // Multithread the warp ['-s_srs', Epsg.get(sourceProjection).toEpsgString()], // Source EPSG @@ -86,38 +57,45 @@ export function gdalBuildVrtWarp( export function gdalBuildCog(targetTiff: URL, sourceVrt: URL, opt: CogifyCreationOptions): GdalCommand { const cfg = { ...Presets[opt.preset], ...opt }; + const tileMatrix = TileMatrixSets.find(cfg.tileMatrix); + if (tileMatrix == null) throw new Error('Unable to find tileMatrix: ' + cfg.tileMatrix); + + const bounds = tileMatrix.tileToSourceBounds(cfg.tile); + const tileExtent = [ + Math.min(bounds.x, bounds.right), + Math.min(bounds.y, bounds.bottom), + Math.max(bounds.x, bounds.right), + Math.max(bounds.y, bounds.bottom), + ]; - const targetOpts = getTargetOptions(cfg); + const targetResolution = tileMatrix.pixelScale(cfg.zoomLevel); return { command: 'gdal_translate', output: targetTiff, args: [ ['-of', 'COG'], - cfg.srcwin ? ['-srcwin', cfg.srcwin[0], cfg.srcwin[1], cfg.srcwin[2], cfg.srcwin[3]] : undefined, - cfg.bigTIFF ? ['-co', `BIGTIFF=${cfg.bigTIFF}`] : ['-co', 'BIGTIFF=IF_NEEDED'], // BigTiff is somewhat slower and most (All?) of the COGS should be well below 4GB ['-co', 'NUM_THREADS=ALL_CPUS'], // Use all CPUS ['--config', 'GDAL_NUM_THREADS', 'all_cpus'], // Also required to NUM_THREADS till gdal 3.7.x + ['-co', 'BIGTIFF=IF_NEEDED'], // BigTiff is somewhat slower and most (All?) of the COGS should be well below 4GB ['-co', 'ADD_ALPHA=YES'], - ['-co', `BLOCKSIZE=${cfg.blockSize}`], /** * GDAL will recompress existing overviews if they exist which will compound * any lossly compression on the overview, so compute new overviews instead */ ['-co', 'OVERVIEWS=IGNORE_EXISTING'], - cfg.overviewCompress ? ['-co', `OVERVIEW_COMPRESS=${cfg.overviewCompress}`] : undefined, - cfg.overviewQuality ? ['-co', `OVERVIEW_QUALITY=${cfg.overviewQuality}`] : undefined, - cfg.warpResampling ? ['-co', `WARP_RESAMPLING=${cfg.warpResampling}`] : undefined, - cfg.overviewResampling ? ['-co', `OVERVIEW_RESAMPLING=${cfg.overviewResampling}`] : undefined, + ['-co', `BLOCKSIZE=${cfg.blockSize}`], + // ['-co', 'RESAMPLING=cubic'], + ['-co', `WARP_RESAMPLING=${cfg.warpResampling}`], + ['-co', `OVERVIEW_RESAMPLING=${cfg.overviewResampling}`], ['-co', `COMPRESS=${cfg.compression}`], cfg.quality ? ['-co', `QUALITY=${cfg.quality}`] : undefined, cfg.maxZError ? ['-co', `MAX_Z_ERROR=${cfg.maxZError}`] : undefined, cfg.maxZErrorOverview ? ['-co', `MAX_Z_ERROR_OVERVIEW=${cfg.maxZErrorOverview}`] : undefined, ['-co', 'SPARSE_OK=YES'], - targetOpts.targetSrs ? ['-co', `TARGET_SRS=${targetOpts.targetSrs}`] : undefined, - targetOpts.extent ? ['-co', `EXTENT=${targetOpts.extent.join(',')},`] : undefined, - targetOpts.targetResolution ? ['-tr', targetOpts.targetResolution, targetOpts.targetResolution] : undefined, - cfg.noReprojecting ? ['-a_srs', `EPSG:${cfg.sourceEpsg}`] : undefined, + ['-co', `TARGET_SRS=${tileMatrix.projection.toEpsgString()}`], + ['-co', `EXTENT=${tileExtent.join(',')},`], + ['-tr', targetResolution, targetResolution], urlToString(sourceVrt), urlToString(targetTiff), ] @@ -143,8 +121,8 @@ export function gdalCreate(targetTiff: URL, color: Rgba, opt: CogifyCreationOpti const tileMatrix = TileMatrixSets.find(cfg.tileMatrix); if (tileMatrix == null) throw new Error('Unable to find tileMatrix: ' + cfg.tileMatrix); - const bounds = tileMatrix.tileToSourceBounds(cfg.tile ?? { x: 0, y: 0, z: 0 }); - const pixelScale = tileMatrix.pixelScale(cfg.zoomLevel ?? 0); + const bounds = tileMatrix.tileToSourceBounds(cfg.tile); + const pixelScale = tileMatrix.pixelScale(cfg.zoomLevel); const size = Math.round(bounds.width / pixelScale); // if the value of 'size' is not a power of 2 @@ -168,3 +146,39 @@ export function gdalCreate(targetTiff: URL, color: Rgba, opt: CogifyCreationOpti .map(String), }; } + +export function gdalBuildTopoRasterCommands(input: URL, output: URL, width: number, height: number): GdalCommand { + const command: GdalCommand = { + output, + command: 'gdal_translate', + args: [ + ['-q'], // Supress non-error output + ['-stats'], // Force stats (re)computation + ['-of', 'COG'], // Output format + ['-srcwin', '0', '0', `${width - DEFAULT_TRIM_PIXEL_RIGHT}`, `${height}`], + + // https://gdal.org/en/latest/drivers/raster/cog.html#creation-options + ['-co', 'BIGTIFF=NO'], + ['-co', 'BLOCKSIZE=512'], + ['-co', 'COMPRESS=WEBP'], + ['-co', 'NUM_THREADS=ALL_CPUS'], // Use all CPUS + ['-co', 'OVERVIEW_COMPRESS=WEBP'], + ['-co', 'OVERVIEWS=IGNORE_EXISTING'], + ['-co', 'OVERVIEW_QUALITY=90'], + ['-co', 'OVERVIEW_RESAMPLING=LANCZOS'], + ['-co', 'QUALITY=100'], + ['-co', 'SPARSE_OK=TRUE'], // Allow for sparse writes + + // https://gdal.org/en/latest/drivers/raster/cog.html#reprojection-related-creation-options + ['-co', 'ADD_ALPHA=YES'], + + urlToString(input), + urlToString(output), + ] + .filter((f) => f != null) + .flat() + .map(String), + }; + + return command; +} diff --git a/packages/cogify/src/cogify/stac.ts b/packages/cogify/src/cogify/stac.ts index 2dfe7e917..cb0f70ce0 100644 --- a/packages/cogify/src/cogify/stac.ts +++ b/packages/cogify/src/cogify/stac.ts @@ -8,11 +8,8 @@ export interface CogifyCreationOptions { /** Preset GDAL config to use */ preset: string; - /** Tile Id to be created */ - tileId: string; - /** Tile to be created */ - tile?: Tile; + tile: Tile; /** Tile matrix to create the tiles against */ tileMatrix: string; @@ -61,36 +58,8 @@ export interface CogifyCreationOptions { */ overviewResampling?: GdalResampling; - /** - * compression method for overview - */ - overviewCompress?: string; - - /** - * JPEG/WEBP quality setting for overviews range from 1 to 100 - */ - overviewQuality?: number; - /** Color with which to replace all transparent COG pixels */ background?: Rgba; - - /** Adds an alpha mask band to the VRT when the source raster have none. */ - addalpha?: boolean; - - /** Stop to reproject the imagery by gdalwarp*/ - noReprojecting?: boolean; - - /** - * External overviews can be created in the BigTIFF format - * - * @default IF_NEEDED - */ - bigTIFF?: 'YES' | 'NO' | 'IF_NEEDED' | 'IF_SAFER'; - - /** - * Selects a subwindow from the source image for copying based on pixel/line location. - */ - srcwin?: number[]; } export type GdalResampling = 'nearest' | 'bilinear' | 'cubic' | 'cubicspline' | 'lanczos' | 'average' | 'mode'; @@ -146,4 +115,4 @@ export function createFileStats(data: string | Buffer): { 'file:size': number; ' // Multihash header for sha256 is 0x12 0x20 'file:checksum': '1220' + createHash('sha256').update(data).digest('hex'), }; -} +} \ No newline at end of file diff --git a/packages/cogify/src/tile.cover.ts b/packages/cogify/src/tile.cover.ts index 7648b48cc..28da4a07a 100644 --- a/packages/cogify/src/tile.cover.ts +++ b/packages/cogify/src/tile.cover.ts @@ -168,7 +168,6 @@ export async function createTileCover(ctx: TileCoverContext): Promise