Update the cogify create cog to support topo raster

Wentao-Kuang · Wentao-Kuang · commit 30904ffde389 · 2025-01-08T14:39:54.000+13:00
diff --git 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, TileId, TileMatrixSet, TileMatrixSets } from '@basemaps/geo';
+import { Projection, ProjectionLoader, 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';
@@ -155,7 +155,7 @@ export const BasemapsCogifyCreateCommand = command({
         const { item, url } = f;
         const cutlineLink = getCutline(item.links);
         const options = item.properties['linz_basemaps:options'];
-        const tileId = TileId.fromTile(options.tile);
+        const tileId = options.tileId;
 
         // Location to where the tiff should be stored
         const tiffPath = new URL(tileId + '.tiff', url);
@@ -268,7 +268,7 @@ export const BasemapsCogifyCreateCommand = command({
       {
         count: toCreate.length,
         created: filtered.length,
-        files: filtered.map((f) => TileId.fromTile(f.item.properties['linz_basemaps:options'].tile)),
+        files: filtered.map((f) => f.item.properties['linz_basemaps:options'].tileId),
       },
       'Cog:Done',
     );
@@ -292,7 +292,7 @@ export interface CogCreationContext {
 async function createCog(ctx: CogCreationContext): Promise<URL> {
   const options = ctx.options;
   await ProjectionLoader.load(options.sourceEpsg);
-  const tileId = TileId.fromTile(options.tile);
+  const tileId = options.tileId;
 
   const logger = ctx.logger?.child({ tileId });
 
@@ -303,13 +303,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);
+  const vrtSourceCommand = gdalBuildVrt(new URL(`${tileId}-source.vrt`, ctx.tempFolder), ctx.sourceFiles, options);
   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) {
+  if (ctx.cutline.path && options.tile) {
     logger?.debug('Cog:Cutline');
     const optimizedCutline = ctx.cutline.optimize(options.tile);
     if (optimizedCutline) {
@@ -321,19 +321,23 @@ async function createCog(ctx: CogCreationContext): Promise<URL> {
     }
   }
 
-  // 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);
+  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;
+  }
 
   if (options.background == null) {
     // Create the COG from the warped vrt without a forced background
-    const cogCreateCommand = gdalBuildCog(new URL(`${tileId}.tiff`, ctx.tempFolder), vrtWarpCommand.output, options);
+    const cogCreateCommand = gdalBuildCog(new URL(`${tileId}.tiff`, ctx.tempFolder), vrtOutput, options);
     await new GdalRunner(cogCreateCommand).run(logger);
     return cogCreateCommand.output;
   }
@@ -345,7 +349,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,
-    vrtWarpCommand.output,
+    vrtOutput,
   ]);
   await new GdalRunner(vrtMergeCommand).run(logger);
 
diff --git a/packages/cogify/src/cogify/cli/cli.cover.ts b/packages/cogify/src/cogify/cli/cli.cover.ts
@@ -144,11 +144,13 @@ export const BasemapsCogifyCoverCommand = command({
     const items = [];
     const tilesByZoom: number[] = [];
     for (const item of res.items) {
-      const tileId = TileId.fromTile(item.properties['linz_basemaps:options'].tile);
+      const tile = item.properties['linz_basemaps:options'].tile;
+      if (tile == null) throw new Error('Tile not found in item');
+      const tileId = TileId.fromTile(tile);
       const itemPath = new URL(`${tileId}.json`, targetPath);
       items.push({ path: itemPath });
       await fsa.write(itemPath, JSON.stringify(item, null, 2));
-      const z = item.properties['linz_basemaps:options'].tile.z;
+      const z = tile.z;
       tilesByZoom[z] = (tilesByZoom[z] ?? 0) + 1;
       ctx.logger?.trace({ path: itemPath }, 'Imagery:Stac:Item:Write');
     }
diff --git a/packages/cogify/src/cogify/gdal.command.ts b/packages/cogify/src/cogify/gdal.command.ts
@@ -8,12 +8,46 @@ import { CogifyCreationOptions } from './stac.js';
 
 const isPowerOfTwo = (x: number): boolean => (x & (x - 1)) === 0;
 
-export function gdalBuildVrt(targetVrt: URL, source: URL[]): GdalCommand {
+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);
+    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 {
   if (source.length === 0) throw new Error('No source files given for :' + targetVrt.href);
   return {
     output: targetVrt,
     command: 'gdalbuildvrt',
-    args: [urlToString(targetVrt), ...source.map(urlToString)],
+    args: [opt && opt.addalpha ? ['-addalpha'] : undefined, urlToString(targetVrt), ...source.map(urlToString)]
+      .filter((f) => f != null)
+      .flat()
+      .map(String),
   };
 }
 
@@ -26,6 +60,7 @@ 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 {
@@ -53,18 +88,8 @@ 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 targetResolution = tileMatrix.pixelScale(cfg.zoomLevel);
+  const targetOpts = getTargetOptions(cfg);
 
   return {
     command: 'gdal_translate',
@@ -73,25 +98,27 @@ export function gdalBuildCog(targetTiff: URL, sourceVrt: URL, opt: CogifyCreatio
       ['-of', 'COG'],
       ['-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
+      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', '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'],
-      ['-co', `BLOCKSIZE=${cfg.blockSize}`],
-      // ['-co', 'RESAMPLING=cubic'],
-      ['-co', `WARP_RESAMPLING=${cfg.warpResampling}`],
-      ['-co', `OVERVIEW_RESAMPLING=${cfg.overviewResampling}`],
+      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', `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'],
-      ['-co', `TARGET_SRS=${tileMatrix.projection.toEpsgString()}`],
-      ['-co', `EXTENT=${tileExtent.join(',')},`],
-      ['-tr', targetResolution, targetResolution],
+      targetOpts.targetSrs ? ['-co', `TARGET_SRS=${targetOpts.targetSrs}`] : undefined,
+      targetOpts.extent ? ['-co', `EXTENT=${targetOpts.extent.join(',')},`] : undefined,
+      targetOpts.targetResolution ? ['-tr', targetOpts.targetResolution, targetOpts.targetResolution] : undefined,
       urlToString(sourceVrt),
       urlToString(targetTiff),
     ]
@@ -117,8 +144,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);
-  const pixelScale = tileMatrix.pixelScale(cfg.zoomLevel);
+  const bounds = tileMatrix.tileToSourceBounds(cfg.tile ?? { x: 0, y: 0, z: 0 });
+  const pixelScale = tileMatrix.pixelScale(cfg.zoomLevel ?? 0);
   const size = Math.round(bounds.width / pixelScale);
 
   // if the value of 'size' is not a power of 2
diff --git a/packages/cogify/src/cogify/stac.ts b/packages/cogify/src/cogify/stac.ts
@@ -8,8 +8,11 @@ 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;
@@ -58,8 +61,36 @@ 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';
diff --git a/packages/cogify/src/tile.cover.ts b/packages/cogify/src/tile.cover.ts
@@ -168,6 +168,7 @@ export async function createTileCover(ctx: TileCoverContext): Promise<TileCoverR
         'linz_basemaps:options': {
           preset: ctx.preset,
           ...Presets[ctx.preset].options,
+          tileId: TileId.fromTile(tile),
           tile,
           tileMatrix: ctx.tileMatrix.identifier,
           sourceEpsg: ctx.imagery.projection,
@@ -225,7 +226,9 @@ export async function createTileCover(ctx: TileCoverContext): Promise<TileCoverR
       temporal: { interval: dateTime.start ? [[dateTime.start, dateTime.end]] : [[cliDate, null]] },
     },
     links: items.map((item) => {
-      const tileId = TileId.fromTile(item.properties['linz_basemaps:options'].tile);
+      const tile = item.properties['linz_basemaps:options'].tile;
+      if (tile == null) throw new Error('Tile missing from item');
+      const tileId = TileId.fromTile(tile);
       return { href: `./${tileId}.json`, rel: 'item', type: 'application/json' };
     }),
   };
