Skip to content

Commit

Permalink
create topo raster cog
Browse files Browse the repository at this point in the history
  • Loading branch information
@Wentao-Kuang
Wentao-Kuang committed Jan 13, 2025
1 parent 07e39b1 commit 87fa033
Show file tree
Hide file tree
Showing 4 changed files with 135 additions and 106 deletions.
95 changes: 71 additions & 24 deletions packages/cogify/src/cogify/cli/cli.cog.ts
View file
Original file line number Diff line number Diff line change
@@ -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';
Expand All @@ -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';
Expand Down Expand Up @@ -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',
Expand Down Expand Up @@ -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);
Expand All @@ -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,
Expand Down Expand Up @@ -268,14 +288,18 @@ 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',
);
},
});

export interface CogCreationContext {
/** TileId for the file name */
tileId: string;
/** COG Creation options */
options: CogifyCreationOptions;
/** Location to store all the temporary files */
Expand All @@ -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;
}
Expand All @@ -292,7 +318,7 @@ export interface CogCreationContext {
async function createCog(ctx: CogCreationContext): Promise<URL> {
const options = ctx.options;
await ProjectionLoader.load(options.sourceEpsg);
const tileId = options.tileId;
const tileId = ctx.tileId;

const logger = ctx.logger?.child({ tileId });

Expand All @@ -303,13 +329,13 @@ async function createCog(ctx: CogCreationContext): Promise<URL> {

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) {
Expand All @@ -321,23 +347,19 @@ async function createCog(ctx: CogCreationContext): Promise<URL> {
}
}

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;
}
Expand All @@ -349,7 +371,7 @@ async function createCog(ctx: CogCreationContext): Promise<URL> {
// 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);

Expand All @@ -359,6 +381,31 @@ async function createCog(ctx: CogCreationContext): Promise<URL> {
return cogCreateCommand.output;
}

/** Create a cog from the creation options */
async function createTopoCog(ctx: CogCreationContext): Promise<URL> {
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
Expand Down
110 changes: 62 additions & 48 deletions packages/cogify/src/cogify/gdal.command.ts
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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),
Expand All @@ -60,14 +30,15 @@ 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 {
output: targetVrt,
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
Expand All @@ -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),
]
Expand All @@ -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
Expand All @@ -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;
}
Loading

0 comments on commit 87fa033

Please sign in to comment.