fix(cli-raster): zstd covering with smaller tiles than expected for single band sources (#3551)

### Motivation

Covering a collection of tiffs with ZSTD was creating more files than
other dataset types.

for some datasets this was a difference in 5000 files

### Modifications

Use the source band information to help determine what target zoom level
a dataset should be tiled at.

### Verification

Unit test + ran some RGB / RGBI tests locally

Author: blacha

packages/cli-raster/src/cogify/covering/__test__/covering.test.ts

@@ -0,0 +1,104 @@
+import assert from 'node:assert';
+import { describe, it } from 'node:test';
+
+import { ImageryBandDataType, ImageryBandType } from '@basemaps/config';
+
+import { findOptimialCoveringZoomOffset } from '../tile.cover.js';
+
+const TileSize = 512;
+
+function rawDataSize(imageryType: ImageryBandDataType): number {
+  switch (imageryType) {
+    case 'uint8':
+      return 1;
+    case 'uint16':
+      return 2;
+    case 'float32':
+      return 4;
+    default:
+      throw new Error(`Unknown imagery type: ${imageryType}`);
+  }
+}
+
+function rawTileSize(sourceBands: ImageryBandType[]): number {
+  let size = 0;
+  for (const b of sourceBands) size += rawDataSize(b.type);
+  return size * TileSize * TileSize;
+}
+
+function tileCount(zoomOffset: number): number {
+  let count = 0;
+  for (let z = 0; z < zoomOffset; z++) {
+    // z:0 - 1x1 - 1 tile
+    // z:1 - 2x2 - 4 tiles
+    // z:3 - 4x4 - 16 tiles
+    // z:3 - 8x8 - 64 tiles
+    count = count + 2 ** z * 2 ** z;
+  }
+  return count;
+}
+
+function estimatedCogSize(zoomOffset: number, sourceBands: ImageryBandType[]): number {
+  const tileCountValue = tileCount(zoomOffset);
+  const tileSizeValue = rawTileSize(sourceBands);
+  return tileCountValue * tileSizeValue;
+}
+
+describe('coveringOffset', () => {
+  const uint8 = { type: 'uint8' } as const;
+  const float32 = { type: 'float32' } as const;
+  const uint16 = { type: 'uint16' } as const;
+
+  // TODO use these estimated sizes to determine if the covering zoom level offset is ok
+  it('should calculate raw size of a 6 level cog correctly', () => {
+    assert.equal(tileCount(1), 1);
+    assert.equal(tileCount(2), 1 + 4);
+    assert.equal(tileCount(3), 1 + 4 + 16);
+    assert.equal(tileCount(4), 1 + 4 + 16 + 64);
+    assert.equal(tileCount(5), 1 + 4 + 16 + 64 + 256);
+
+    assert.equal(estimatedCogSize(1, [uint8]), 1 * 512 * 512 * 1);
+    assert.equal(estimatedCogSize(1, [uint16]), 1 * 512 * 512 * 2);
+    assert.equal(estimatedCogSize(1, [uint16, uint16]), 1 * 512 * 512 * 2 * 2);
+    assert.equal(estimatedCogSize(1, [float32]), 1 * 512 * 512 * 4);
+
+    rawTileSize([uint16, uint16, uint16, uint16, uint16]);
+  });
+
+  it('should cover singleband uint8 or uint16 correctly', () => {
+    for (const preset of ['lerc_1mm', 'lerc_10mm', 'zstd_17'] as const) {
+      assert.equal(7, findOptimialCoveringZoomOffset(preset, [uint8]));
+      assert.equal(7, findOptimialCoveringZoomOffset(preset, [uint16]));
+    }
+  });
+
+  it('should cover a webp rgb(a) correctly', () => {
+    assert.equal(7, findOptimialCoveringZoomOffset('webp', [uint8, uint8, uint8, uint8]));
+    assert.equal(7, findOptimialCoveringZoomOffset('webp', [uint8, uint8, uint8]));
+  });
+
+  it('should cover a float32 single band correctly', () => {
+    assert.equal(7, findOptimialCoveringZoomOffset('lerc_1mm', [float32]));
+    assert.equal(7, findOptimialCoveringZoomOffset('lerc_10mm', [float32]));
+  });
+
+  it('should reduce the size with zstd compression', () => {
+    assert.equal(6, findOptimialCoveringZoomOffset('zstd_17', [float32]));
+    assert.equal(6, findOptimialCoveringZoomOffset('zstd_17', [uint8, uint8, uint8]));
+    assert.equal(6, findOptimialCoveringZoomOffset('zstd_17', [uint8, uint8, uint8, uint8]));
+
+    // RGBI+Alpha
+    assert.equal(6, findOptimialCoveringZoomOffset('zstd_17', [uint8, uint8, uint8, uint8, uint8]));
+    assert.equal(6, findOptimialCoveringZoomOffset('zstd_17', [uint16, uint16, uint16, uint16, uint16]));
+  });
+
+  it('should default to level 6 for lzw compression', () => {
+    assert.equal(6, findOptimialCoveringZoomOffset('lzw', [float32]));
+    assert.equal(6, findOptimialCoveringZoomOffset('lzw', [uint8, uint8, uint8]));
+    assert.equal(6, findOptimialCoveringZoomOffset('lzw', [uint8, uint8, uint8, uint8]));
+
+    // RGBI+Alpha - these have not been tested and are assumed to be ok
+    assert.equal(6, findOptimialCoveringZoomOffset('lzw', [uint8, uint8, uint8, uint8, uint8]));
+    assert.equal(6, findOptimialCoveringZoomOffset('lzw', [uint16, uint16, uint16, uint16, uint16]));
+  });
+});

packages/cli-raster/src/cogify/covering/tile.cover.ts

@@ -1,4 +1,4 @@
-import { Rgba } from '@basemaps/config';
+import { ImageryBandType, Rgba } from '@basemaps/config';
 import { ConfigImageryTiff } from '@basemaps/config-loader';
 import { BoundingBox, Bounds, EpsgCode, Projection, ProjectionLoader, TileId, TileMatrixSet } from '@basemaps/geo';
 import { fsa, LogType, urlToString } from '@basemaps/shared';
@@ -75,17 +75,44 @@ function getTargetBaseZoom(tileMatrix: TileMatrixSet, resolution: number, target
   return Projection.getTiffResZoom(tileMatrix, resolution) + targetZoomOffset;
 }
 
-// The base zoom is 256x256 pixels at its resolution, we are trying to find a image that is <32k pixels wide/high
-// zooming out 7 levels converts a 256x256 image into 32k x 32k image
-// 256 * 2 ** 7 = 32,768 - 256x256 tile
-// 512 * 2 ** 6 = 32,768 - 512x512 tile
-// This math only works for highly compressed RGB imagery, for multispectrial imagery small tiles need to be made
 export const TargetZoomOffsetDefault = 7;
-// ZSTD files are generally larger than webp or LERC
-export const TargetZoomOffset: Partial<Record<PresetName, number>> = {
-  zstd_17: 6,
-  lzw: 6,
-};
+
+/**
+ * We are trying to find a zoom level that will create COGS that are optimally sized for the preset,
+ * We are looking for cogs that are around 1GB in size on average and keeping under 4GB when maximally filled
+ *
+ * There are lots of different presets and source imagery bands which affect the final size of the COGs,
+ * for examples:
+ * - RGB(A) `uint8` imagery compressed with webp this preset is approximately 7 levels from the base zoom level.
+ * - Elevation `float32` compressed with zstd is less effective the offset is reduced to 6 levels
+ *
+ * @param preset target compression preset
+ * @param targetBaseZoom the base zoom level to use
+ * @param sourceBands bands present in the source dataset
+ *
+ * @returns a zoom level offset to use for the covering
+ */
+export function findOptimialCoveringZoomOffset(preset: PresetName, sourceBands: ImageryBandType[]): number {
+  // Webp is very efficient at compressing RGB(A) imagery so we can keep the default offset
+  if (preset === 'webp') return TargetZoomOffsetDefault;
+  // LZW is not very effective at compressing most imagery so we reduce the offset by one
+  if (preset === 'lzw') return TargetZoomOffsetDefault - 1;
+
+  const sourceBandText = sourceBands.map((m) => m.type).join(',');
+  // Single band float32 imagery tends to be elevation data which compresses well with lossy lerc and less so with lossless zstd
+  if (sourceBandText === 'float32') {
+    if (preset === 'lerc_1mm' || preset === 'lerc_10mm') return TargetZoomOffsetDefault;
+    // all other compressors are lossless
+    return TargetZoomOffsetDefault - 1;
+  }
+
+  // Single band uint8 or uint16 compresses generally very effectively with most presets
+  if (sourceBandText === 'uint8') return TargetZoomOffsetDefault;
+  if (sourceBandText === 'uint16') return TargetZoomOffsetDefault;
+
+  // Unknown preset so assume one level offset is ok
+  return TargetZoomOffsetDefault - 1;
+}
 
 export async function createTileCover(ctx: TileCoverContext): Promise<TileCoverResult> {
   // Ensure we have the projection loaded for the source imagery
@@ -94,10 +121,10 @@ export async function createTileCover(ctx: TileCoverContext): Promise<TileCoverR
 
   // Find the zoom level that is at least as good as the source imagery
   const targetBaseZoom = getTargetBaseZoom(ctx.tileMatrix, ctx.imagery.gsd, ctx.targetZoomOffset);
-
-  const targetZoomOffset = TargetZoomOffset[ctx.preset] ?? 7;
-
-  const optimalCoveringZoom = Math.max(1, targetBaseZoom - targetZoomOffset); // z12 from z19
+  const optimalCoveringZoom = Math.max(
+    1,
+    targetBaseZoom - findOptimialCoveringZoomOffset(ctx.preset, ctx.imagery.bands),
+  );
   ctx.logger?.debug({ targetBaseZoom, cogOverZoom: optimalCoveringZoom }, 'Imagery:ZoomLevel');
 
   const sourceBounds = projectPolygon(