From 753da7e5f5c6634e989cd3497bb0cef2b6eecde2 Mon Sep 17 00:00:00 2001 From: Blayne Chard Date: Tue, 15 Oct 2024 16:50:21 +1300 Subject: [PATCH 01/12] wip: retile datasets --- .../tileindex-validate/tileindex.validate.ts | 91 ++++++++++++++----- 1 file changed, 66 insertions(+), 25 deletions(-) diff --git a/src/commands/tileindex-validate/tileindex.validate.ts b/src/commands/tileindex-validate/tileindex.validate.ts index f97c5e5b..7e929f59 100644 --- a/src/commands/tileindex-validate/tileindex.validate.ts +++ b/src/commands/tileindex-validate/tileindex.validate.ts @@ -1,6 +1,9 @@ -import { Bounds, Projection } from '@basemaps/geo'; +import assert from 'node:assert'; + +import { Bounds, Point, Projection } from '@basemaps/geo'; import { fsa } from '@chunkd/fs'; import { Size, Tiff, TiffTag } from '@cogeotiff/core'; +import { BBox } from '@linzjs/geojson'; import { boolean, command, flag, number, option, optional, restPositionals, string, Type } from 'cmd-ts'; import { CliInfo } from '../../cli.info.js'; @@ -242,8 +245,8 @@ export const commandTileIndexValidate = command({ } return Projection.get(epsg).boundsToGeoJsonFeature(Bounds.fromBbox(loc.bbox), { source: loc.source, - tileName: loc.tileName, - isDuplicate: (outputs.get(loc.tileName)?.length ?? 1) > 1, + tileName: loc.tileNames.join(', '), + isDuplicate: true, // (outputs.get(loc.tileNames)?.length ?? 1) > 1, }); }), }); @@ -251,14 +254,12 @@ export const commandTileIndexValidate = command({ await fsa.write('/tmp/tile-index-validate/output.geojson', { type: 'FeatureCollection', - features: [...outputs.values()].map((locs) => { - const firstLoc = locs[0]; - if (firstLoc == null) throw new Error('Unable to extract tiff locations from: ' + args.location.join(', ')); - const mapTileIndex = MapSheet.getMapTileIndex(firstLoc.tileName); - if (mapTileIndex == null) throw new Error('Failed to extract tile information from: ' + firstLoc.tileName); + features: [...outputs.keys()].map((key) => { + const mapTileIndex = MapSheet.getMapTileIndex(key); + if (mapTileIndex == null) throw new Error('Failed to extract tile information from: ' + key); return Projection.get(2193).boundsToGeoJsonFeature(Bounds.fromBbox(mapTileIndex.bbox), { - source: locs.map((l) => l.source), - tileName: firstLoc.tileName, + source: outputs.get(key)?.map((l) => l.source), + tileName: key, }); }), }); @@ -266,8 +267,9 @@ export const commandTileIndexValidate = command({ await fsa.write( '/tmp/tile-index-validate/file-list.json', - [...outputs.values()].map((locs) => { - return { output: locs[0]?.tileName, input: locs.map((l) => l.source) }; + [...outputs.keys()].map((key) => { + const locs = outputs.get(key); + return { output: key, input: locs?.map((l) => l.source) }; }), ); logger.info({ path: '/tmp/tile-index-validate/file-list.json', count: outputs.size }, 'Write:FileList'); @@ -333,9 +335,11 @@ function validateConsistentBands(locs: TiffLocation[]): string[] { export function groupByTileName(tiffs: TiffLocation[]): Map { const duplicates: Map = new Map(); for (const loc of tiffs) { - const uris = duplicates.get(loc.tileName) ?? []; - uris.push(loc); - duplicates.set(loc.tileName, uris); + for (const sheetCode of loc.tileNames) { + const uris = duplicates.get(sheetCode) ?? []; + uris.push(loc); + duplicates.set(sheetCode, uris); + } } return duplicates; } @@ -348,7 +352,7 @@ export interface TiffLocation { /** EPSG code of the tiff if found */ epsg?: number | null; /** Output tile name */ - tileName: string; + tileNames: string[]; /** * List of bands inside the tiff in the format `uint8` `uint16` * @@ -399,6 +403,14 @@ export async function extractTiffLocations( // Tilename from center const tileName = getTileName(x, y, gridSize); + const [ulX, ulY] = targetProjection.fromWgs84(sourceProjection.toWgs84([bbox[0], bbox[3]])); + const [lrX, lrY] = targetProjection.fromWgs84(sourceProjection.toWgs84([bbox[2], bbox[1]])); + + console.log('topLeft', MapSheet.sheetCode(ulX!, ulY!)); + + const covering = [...iterateMapSheets([ulX, ulY, lrX, lrY], gridSize)]; + assert.ok(covering.includes(tileName)); + // if (shouldValidate) { // // Is the tiff bounding box the same as the map sheet bounding box! // // Also need to allow for ~1.5cm of error between bounding boxes. @@ -407,7 +419,7 @@ export async function extractTiffLocations( return { bbox, source: tiff.source.url.href, - tileName, + tileNames: covering, epsg: tiff.images[0]?.epsg, bands: await extractBandInformation(tiff), }; @@ -434,10 +446,13 @@ export function getSize(extent: [number, number, number, number]): Size { } export function validateTiffAlignment(tiff: TiffLocation, allowedError = 0.015): boolean { - const mapTileIndex = MapSheet.getMapTileIndex(tiff.tileName); + const tileName = tiff.tileNames[0]; + // FIXME + if (tileName == null) return false; + const mapTileIndex = MapSheet.getMapTileIndex(tileName); if (mapTileIndex == null) { logger.error( - { reason: `Failed to extract bounding box from: ${tiff.tileName}`, source: tiff.source }, + { reason: `Failed to extract bounding box from: ${tileName}`, source: tiff.source }, 'TileInvalid:Validation:Failed', ); return false; @@ -478,12 +493,7 @@ export function validateTiffAlignment(tiff: TiffLocation, allowedError = 0.015): } export function getTileName(x: number, y: number, gridSize: GridSize): string { - const offsetX = Math.round(Math.floor((x - MapSheet.origin.x) / MapSheet.width)); - const offsetY = Math.round(Math.floor((MapSheet.origin.y - y) / MapSheet.height)); - - // Build name - const letters = Object.keys(SheetRanges)[offsetY]; - const sheetCode = `${letters}${`${offsetX}`.padStart(2, '0')}`; + const sheetCode = MapSheet.sheetCode(x, y); // TODO: re-enable this check when validation logic // if (!MapSheet.isKnown(sheetCode)) throw new Error('Map sheet outside known range: ' + sheetCode); @@ -496,6 +506,8 @@ export function getTileName(x: number, y: number, gridSize: GridSize): string { const nbDigits = gridSize === 500 ? 3 : 2; + const offsetX = Math.round(Math.floor((x - MapSheet.origin.x) / MapSheet.width)); + const offsetY = Math.round(Math.floor((MapSheet.origin.y - y) / MapSheet.height)); const maxY = MapSheet.origin.y - offsetY * MapSheet.height; const minX = MapSheet.origin.x + offsetX * MapSheet.width; const tileX = Math.round(Math.floor((x - minX) / tileWidth + 1)); @@ -522,3 +534,32 @@ export async function validate8BitsTiff(tiff: Tiff): Promise { throw new Error(`${tiff.source.url.href} is not a 8 bits TIFF`); } } + +export function* iterateMapSheets(bounds: BBox, gridSize: GridSize): Iterator { + const minX = Math.min(bounds[0], bounds[2]); + const maxX = Math.max(bounds[0], bounds[2]); + const minY = Math.min(bounds[1], bounds[3]); + const maxY = Math.max(bounds[1], bounds[3]); + + // const minOffsetX = Math.round(Math.floor((minX - MapSheet.origin.x) / MapSheet.width)); + // const minOffsetY = Math.round(Math.floor((MapSheet.origin.y - maxY) / MapSheet.height)); + + // const maxOffsetX = Math.round(Math.floor((maxX - MapSheet.origin.x) / MapSheet.width)); + // const maxOffsetY = Math.round(Math.floor((MapSheet.origin.y - minY) / MapSheet.height)); + // console.log({ minOffsetX, minOffsetY }, { maxOffsetX, maxOffsetY }); + + // const offsetX = Math.round(Math.floor((minX - MapSheet.origin.x) / MapSheet.width)); + // const offsetY = Math.round(Math.floor((MapSheet.origin.y - minY) / MapSheet.height)); + // console.log(bounds); + + const tilesPerMapSheet = Math.floor(MapSheet.gridSizeMax / gridSize); + console.log({ minX, minY, maxX, maxY }); + + const tileWidth = Math.floor(MapSheet.width / tilesPerMapSheet); + const tileHeight = Math.floor(MapSheet.height / tilesPerMapSheet); + for (let x = minX; x <= maxX + tileWidth - 1; x += tileWidth) { + for (let y = maxY; y >= minY - tileHeight - 1; y -= tileHeight) { + yield getTileName(x, y, gridSize); + } + } +} From e493fcd8a307c1ea1d6818964261636e8036f743 Mon Sep 17 00:00:00 2001 From: Blayne Chard Date: Tue, 15 Oct 2024 16:51:53 +1300 Subject: [PATCH 02/12] wip: refactor build --- .../__test__/tileindex.validate.test.ts | 6 +++--- .../tileindex-validate/tileindex.validate.ts | 21 ++++++++++++------- 2 files changed, 17 insertions(+), 10 deletions(-) diff --git a/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts b/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts index f182c88d..4a060dfe 100644 --- a/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts +++ b/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts @@ -89,8 +89,8 @@ describe('tiffLocation', () => { TiffAy29.images[0].origin[0] = 1684000; TiffAy29.images[0].origin[1] = 6018000; const location = await extractTiffLocations([TiffAs21, TiffAy29], 1000); - assert.equal(location[0]?.tileName, 'AS21_1000_0101'); - assert.equal(location[1]?.tileName, 'AY29_1000_0101'); + assert.equal(location[0]?.tileNames[0], 'AS21_1000_0101'); + assert.equal(location[1]?.tileNames[0], 'AY29_1000_0101'); }); it('should find duplicates', async () => { @@ -118,7 +118,7 @@ describe('tiffLocation', () => { TiffAy29.images[0].origin[0] = 19128043.69337794; TiffAy29.images[0].origin[1] = -4032710.6009459053; const location = await extractTiffLocations([TiffAy29], 1000); - assert.equal(location[0]?.tileName, 'AS21_1000_0101'); + assert.equal(location[0]?.tileNames[0], 'AS21_1000_0101'); }); it('should fail if one location is not extracted', async () => { diff --git a/src/commands/tileindex-validate/tileindex.validate.ts b/src/commands/tileindex-validate/tileindex.validate.ts index 7e929f59..ce56b113 100644 --- a/src/commands/tileindex-validate/tileindex.validate.ts +++ b/src/commands/tileindex-validate/tileindex.validate.ts @@ -1,6 +1,6 @@ import assert from 'node:assert'; -import { Bounds, Point, Projection } from '@basemaps/geo'; +import { Bounds, Projection } from '@basemaps/geo'; import { fsa } from '@chunkd/fs'; import { Size, Tiff, TiffTag } from '@cogeotiff/core'; import { BBox } from '@linzjs/geojson'; @@ -12,7 +12,7 @@ import { isArgo } from '../../utils/argo.js'; import { extractBandInformation } from '../../utils/band.js'; import { FileFilter, getFiles } from '../../utils/chunk.js'; import { findBoundingBox } from '../../utils/geotiff.js'; -import { GridSize, GridSizes, MapSheet, MapSheetTileGridSize, SheetRanges } from '../../utils/mapsheet.js'; +import { GridSize, GridSizes, MapSheet, MapSheetTileGridSize } from '../../utils/mapsheet.js'; import { config, createTiff, forceOutput, registerCli, verbose } from '../common.js'; import { CommandListArgs } from '../list/list.js'; @@ -278,15 +278,16 @@ export const commandTileIndexValidate = command({ let retileNeeded = false; for (const val of outputs.values()) { if (val.length < 2) continue; + // FIXME if (args.retile) { const bandType = validateConsistentBands(val); logger.info( - { tileName: val[0]?.tileName, uris: val.map((v) => v.source), bands: bandType }, + { tileName: val[0]?.tileNames[0], uris: val.map((v) => v.source), bands: bandType }, 'TileIndex:Retile', ); } else { retileNeeded = true; - logger.error({ tileName: val[0]?.tileName, uris: val.map((v) => v.source) }, 'TileIndex:Duplicate'); + logger.error({ tileName: val[0]?.tileNames[0], uris: val.map((v) => v.source) }, 'TileIndex:Duplicate'); } } @@ -406,9 +407,15 @@ export async function extractTiffLocations( const [ulX, ulY] = targetProjection.fromWgs84(sourceProjection.toWgs84([bbox[0], bbox[3]])); const [lrX, lrY] = targetProjection.fromWgs84(sourceProjection.toWgs84([bbox[2], bbox[1]])); - console.log('topLeft', MapSheet.sheetCode(ulX!, ulY!)); + if (ulX == null || ulY == null || lrX == null || lrY == null) { + logger.error( + { reason: 'Failed to reproject point', source: tiff.source }, + 'Reprojection:ExtracTiffLocations:Failed', + ); + return null; + } - const covering = [...iterateMapSheets([ulX, ulY, lrX, lrY], gridSize)]; + const covering = [...iterateMapSheets([ulX, ulY, lrX, lrY], gridSize)] as string[]; assert.ok(covering.includes(tileName)); // if (shouldValidate) { @@ -535,7 +542,7 @@ export async function validate8BitsTiff(tiff: Tiff): Promise { } } -export function* iterateMapSheets(bounds: BBox, gridSize: GridSize): Iterator { +export function* iterateMapSheets(bounds: BBox, gridSize: GridSize): Generator { const minX = Math.min(bounds[0], bounds[2]); const maxX = Math.max(bounds[0], bounds[2]); const minY = Math.min(bounds[1], bounds[3]); From ed5c4e0a23a0857934e85c1ea9e1eb0d6aa4a857 Mon Sep 17 00:00:00 2001 From: Blayne Chard Date: Tue, 15 Oct 2024 16:53:39 +1300 Subject: [PATCH 03/12] refactor: broken test --- .../__test__/tileindex.validate.test.ts | 58 +++++++++---------- 1 file changed, 29 insertions(+), 29 deletions(-) diff --git a/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts b/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts index 4a060dfe..7ad6c8b5 100644 --- a/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts +++ b/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts @@ -150,38 +150,38 @@ describe('validate', () => { sourceEpsg: undefined, }; - it('should fail if duplicate tiles are detected', async (t) => { - // Input source/a/AS21_1000_0101.tiff source/b/AS21_1000_0101.tiff - const stub = t.mock.method(TiffLoader, 'load', () => - Promise.resolve([FakeCogTiff.fromTileName('AS21_1000_0101'), FakeCogTiff.fromTileName('AS21_1000_0101')]), - ); + // it('should fail if duplicate tiles are detected', async (t) => { + // // Input source/a/AS21_1000_0101.tiff source/b/AS21_1000_0101.tiff + // const stub = t.mock.method(TiffLoader, 'load', () => + // Promise.resolve([FakeCogTiff.fromTileName('AS21_1000_0101'), FakeCogTiff.fromTileName('AS21_1000_0101')]), + // ); - try { - await commandTileIndexValidate.handler({ - ...baseArguments, - location: ['s3://test'], - retile: false, - validate: true, - scale: 1000, - forceOutput: true, - }); - assert.fail('Should throw exception'); - } catch (e) { - assert.equal(String(e), 'Error: Duplicate files found, see output.geojson'); - } + // try { + // await commandTileIndexValidate.handler({ + // ...baseArguments, + // location: ['s3://test'], + // retile: false, + // validate: true, + // scale: 1000, + // forceOutput: true, + // }); + // assert.fail('Should throw exception'); + // } catch (e) { + // assert.equal(String(e), 'Error: Duplicate files found, see output.geojson'); + // } - assert.equal(stub.mock.callCount(), 1); - assert.deepEqual(stub.mock.calls[0]?.arguments[0], ['s3://test']); + // assert.equal(stub.mock.callCount(), 1); + // assert.deepEqual(stub.mock.calls[0]?.arguments[0], ['s3://test']); - const outputFileList: FeatureCollection = await fsa.readJson('/tmp/tile-index-validate/output.geojson'); - assert.equal(outputFileList.features.length, 1); - const firstFeature = outputFileList.features[0]; - assert.equal(firstFeature?.properties?.['tileName'], 'AS21_1000_0101'); - assert.deepEqual(firstFeature?.properties?.['source'], [ - 's3://path/AS21_1000_0101.tiff', - 's3://path/AS21_1000_0101.tiff', - ]); - }); + // const outputFileList: FeatureCollection = await fsa.readJson('/tmp/tile-index-validate/output.geojson'); + // assert.equal(outputFileList.features.length, 1); + // const firstFeature = outputFileList.features[0]; + // assert.equal(firstFeature?.properties?.['tileName'], 'AS21_1000_0101'); + // assert.deepEqual(firstFeature?.properties?.['source'], [ + // 's3://path/AS21_1000_0101.tiff', + // 's3://path/AS21_1000_0101.tiff', + // ]); + // }); it('should not fail if duplicate tiles are detected but --retile is used', async (t) => { // Input source/a/AS21_1000_0101.tiff source/b/AS21_1000_0101.tiff From df37e4ac7988b4ef161c979d0085a3a03edcac36 Mon Sep 17 00:00:00 2001 From: Blayne Chard Date: Tue, 15 Oct 2024 16:54:57 +1300 Subject: [PATCH 04/12] refactor: disable more tests --- .../__test__/tileindex.validate.test.ts | 35 +++++++++---------- 1 file changed, 17 insertions(+), 18 deletions(-) diff --git a/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts b/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts index 7ad6c8b5..b7767e80 100644 --- a/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts +++ b/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts @@ -3,7 +3,6 @@ import { before, beforeEach, describe, it } from 'node:test'; import { fsa } from '@chunkd/fs'; import { FsMemory } from '@chunkd/source-memory'; -import { FeatureCollection } from 'geojson'; import { MapSheetData } from '../../../utils/__test__/mapsheet.data.js'; import { GridSize, MapSheet } from '../../../utils/mapsheet.js'; @@ -183,24 +182,24 @@ describe('validate', () => { // ]); // }); - it('should not fail if duplicate tiles are detected but --retile is used', async (t) => { - // Input source/a/AS21_1000_0101.tiff source/b/AS21_1000_0101.tiff - t.mock.method(TiffLoader, 'load', () => - Promise.resolve([FakeCogTiff.fromTileName('AS21_1000_0101'), FakeCogTiff.fromTileName('AS21_1000_0101')]), - ); + // it('should not fail if duplicate tiles are detected but --retile is used', async (t) => { + // // Input source/a/AS21_1000_0101.tiff source/b/AS21_1000_0101.tiff + // t.mock.method(TiffLoader, 'load', () => + // Promise.resolve([FakeCogTiff.fromTileName('AS21_1000_0101'), FakeCogTiff.fromTileName('AS21_1000_0101')]), + // ); - await commandTileIndexValidate.handler({ - ...baseArguments, - location: ['s3://test'], - retile: true, - scale: 1000, - forceOutput: true, - }); - const outputFileList = await fsa.readJson('/tmp/tile-index-validate/file-list.json'); - assert.deepEqual(outputFileList, [ - { output: 'AS21_1000_0101', input: ['s3://path/AS21_1000_0101.tiff', 's3://path/AS21_1000_0101.tiff'] }, - ]); - }); + // await commandTileIndexValidate.handler({ + // ...baseArguments, + // location: ['s3://test'], + // retile: true, + // scale: 1000, + // forceOutput: true, + // }); + // const outputFileList = await fsa.readJson('/tmp/tile-index-validate/file-list.json'); + // assert.deepEqual(outputFileList, [ + // { output: 'AS21_1000_0101', input: ['s3://path/AS21_1000_0101.tiff', 's3://path/AS21_1000_0101.tiff'] }, + // ]); + // }); for (const offset of [0.05, -0.05]) { it(`should fail if input tiff origin X is offset by ${offset}m`, async (t) => { From a7a28c7c949fd069a03f79835701f5c1f65009ec Mon Sep 17 00:00:00 2001 From: Tobias Schmidt <13055656+schmidtnz@users.noreply.github.com> Date: Tue, 10 Dec 2024 11:32:24 +1300 Subject: [PATCH 05/12] fix: iterateMapSheets() returns too many tiles --- src/commands/tileindex-validate/tileindex.validate.ts | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/commands/tileindex-validate/tileindex.validate.ts b/src/commands/tileindex-validate/tileindex.validate.ts index ce56b113..0b6fa085 100644 --- a/src/commands/tileindex-validate/tileindex.validate.ts +++ b/src/commands/tileindex-validate/tileindex.validate.ts @@ -564,8 +564,8 @@ export function* iterateMapSheets(bounds: BBox, gridSize: GridSize): Generator= minY - tileHeight - 1; y -= tileHeight) { + for (let x = minX; x < maxX; x += tileWidth) { + for (let y = maxY; y > minY; y -= tileHeight) { yield getTileName(x, y, gridSize); } } From bddd9e1b05d9243bf96218182a5b5e4ce4564940 Mon Sep 17 00:00:00 2001 From: Tobias Schmidt <13055656+schmidtnz@users.noreply.github.com> Date: Tue, 10 Dec 2024 11:33:21 +1300 Subject: [PATCH 06/12] test: re-enabled tests --- .../__test__/tileindex.validate.test.ts | 93 ++++++++++--------- 1 file changed, 47 insertions(+), 46 deletions(-) diff --git a/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts b/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts index b7767e80..4a060dfe 100644 --- a/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts +++ b/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts @@ -3,6 +3,7 @@ import { before, beforeEach, describe, it } from 'node:test'; import { fsa } from '@chunkd/fs'; import { FsMemory } from '@chunkd/source-memory'; +import { FeatureCollection } from 'geojson'; import { MapSheetData } from '../../../utils/__test__/mapsheet.data.js'; import { GridSize, MapSheet } from '../../../utils/mapsheet.js'; @@ -149,57 +150,57 @@ describe('validate', () => { sourceEpsg: undefined, }; - // it('should fail if duplicate tiles are detected', async (t) => { - // // Input source/a/AS21_1000_0101.tiff source/b/AS21_1000_0101.tiff - // const stub = t.mock.method(TiffLoader, 'load', () => - // Promise.resolve([FakeCogTiff.fromTileName('AS21_1000_0101'), FakeCogTiff.fromTileName('AS21_1000_0101')]), - // ); + it('should fail if duplicate tiles are detected', async (t) => { + // Input source/a/AS21_1000_0101.tiff source/b/AS21_1000_0101.tiff + const stub = t.mock.method(TiffLoader, 'load', () => + Promise.resolve([FakeCogTiff.fromTileName('AS21_1000_0101'), FakeCogTiff.fromTileName('AS21_1000_0101')]), + ); - // try { - // await commandTileIndexValidate.handler({ - // ...baseArguments, - // location: ['s3://test'], - // retile: false, - // validate: true, - // scale: 1000, - // forceOutput: true, - // }); - // assert.fail('Should throw exception'); - // } catch (e) { - // assert.equal(String(e), 'Error: Duplicate files found, see output.geojson'); - // } + try { + await commandTileIndexValidate.handler({ + ...baseArguments, + location: ['s3://test'], + retile: false, + validate: true, + scale: 1000, + forceOutput: true, + }); + assert.fail('Should throw exception'); + } catch (e) { + assert.equal(String(e), 'Error: Duplicate files found, see output.geojson'); + } - // assert.equal(stub.mock.callCount(), 1); - // assert.deepEqual(stub.mock.calls[0]?.arguments[0], ['s3://test']); + assert.equal(stub.mock.callCount(), 1); + assert.deepEqual(stub.mock.calls[0]?.arguments[0], ['s3://test']); - // const outputFileList: FeatureCollection = await fsa.readJson('/tmp/tile-index-validate/output.geojson'); - // assert.equal(outputFileList.features.length, 1); - // const firstFeature = outputFileList.features[0]; - // assert.equal(firstFeature?.properties?.['tileName'], 'AS21_1000_0101'); - // assert.deepEqual(firstFeature?.properties?.['source'], [ - // 's3://path/AS21_1000_0101.tiff', - // 's3://path/AS21_1000_0101.tiff', - // ]); - // }); + const outputFileList: FeatureCollection = await fsa.readJson('/tmp/tile-index-validate/output.geojson'); + assert.equal(outputFileList.features.length, 1); + const firstFeature = outputFileList.features[0]; + assert.equal(firstFeature?.properties?.['tileName'], 'AS21_1000_0101'); + assert.deepEqual(firstFeature?.properties?.['source'], [ + 's3://path/AS21_1000_0101.tiff', + 's3://path/AS21_1000_0101.tiff', + ]); + }); - // it('should not fail if duplicate tiles are detected but --retile is used', async (t) => { - // // Input source/a/AS21_1000_0101.tiff source/b/AS21_1000_0101.tiff - // t.mock.method(TiffLoader, 'load', () => - // Promise.resolve([FakeCogTiff.fromTileName('AS21_1000_0101'), FakeCogTiff.fromTileName('AS21_1000_0101')]), - // ); + it('should not fail if duplicate tiles are detected but --retile is used', async (t) => { + // Input source/a/AS21_1000_0101.tiff source/b/AS21_1000_0101.tiff + t.mock.method(TiffLoader, 'load', () => + Promise.resolve([FakeCogTiff.fromTileName('AS21_1000_0101'), FakeCogTiff.fromTileName('AS21_1000_0101')]), + ); - // await commandTileIndexValidate.handler({ - // ...baseArguments, - // location: ['s3://test'], - // retile: true, - // scale: 1000, - // forceOutput: true, - // }); - // const outputFileList = await fsa.readJson('/tmp/tile-index-validate/file-list.json'); - // assert.deepEqual(outputFileList, [ - // { output: 'AS21_1000_0101', input: ['s3://path/AS21_1000_0101.tiff', 's3://path/AS21_1000_0101.tiff'] }, - // ]); - // }); + await commandTileIndexValidate.handler({ + ...baseArguments, + location: ['s3://test'], + retile: true, + scale: 1000, + forceOutput: true, + }); + const outputFileList = await fsa.readJson('/tmp/tile-index-validate/file-list.json'); + assert.deepEqual(outputFileList, [ + { output: 'AS21_1000_0101', input: ['s3://path/AS21_1000_0101.tiff', 's3://path/AS21_1000_0101.tiff'] }, + ]); + }); for (const offset of [0.05, -0.05]) { it(`should fail if input tiff origin X is offset by ${offset}m`, async (t) => { From 9657405245261203716e432a1a919180731ed22a Mon Sep 17 00:00:00 2001 From: Tobias Schmidt <13055656+schmidtnz@users.noreply.github.com> Date: Tue, 10 Dec 2024 13:59:32 +1300 Subject: [PATCH 07/12] fix: apply rounding to reprojected tile bounds --- .../tileindex-validate/tileindex.validate.ts | 20 ++++--------------- 1 file changed, 4 insertions(+), 16 deletions(-) diff --git a/src/commands/tileindex-validate/tileindex.validate.ts b/src/commands/tileindex-validate/tileindex.validate.ts index 0b6fa085..f3309a84 100644 --- a/src/commands/tileindex-validate/tileindex.validate.ts +++ b/src/commands/tileindex-validate/tileindex.validate.ts @@ -392,7 +392,7 @@ export async function extractTiffLocations( const targetProjection = Projection.get(2193); const sourceProjection = Projection.get(sourceEpsg); - const [x, y] = targetProjection.fromWgs84(sourceProjection.toWgs84([centerX, centerY])); + const [x, y] = targetProjection.fromWgs84(sourceProjection.toWgs84([centerX, centerY])).map(Math.round); if (x == null || y == null) { logger.error( { reason: 'Failed to reproject point', source: tiff.source }, @@ -404,8 +404,8 @@ export async function extractTiffLocations( // Tilename from center const tileName = getTileName(x, y, gridSize); - const [ulX, ulY] = targetProjection.fromWgs84(sourceProjection.toWgs84([bbox[0], bbox[3]])); - const [lrX, lrY] = targetProjection.fromWgs84(sourceProjection.toWgs84([bbox[2], bbox[1]])); + const [ulX, ulY] = targetProjection.fromWgs84(sourceProjection.toWgs84([bbox[0], bbox[3]])).map(Math.round); + const [lrX, lrY] = targetProjection.fromWgs84(sourceProjection.toWgs84([bbox[2], bbox[1]])).map(Math.round); if (ulX == null || ulY == null || lrX == null || lrY == null) { logger.error( @@ -502,7 +502,7 @@ export function validateTiffAlignment(tiff: TiffLocation, allowedError = 0.015): export function getTileName(x: number, y: number, gridSize: GridSize): string { const sheetCode = MapSheet.sheetCode(x, y); // TODO: re-enable this check when validation logic - // if (!MapSheet.isKnown(sheetCode)) throw new Error('Map sheet outside known range: ' + sheetCode); + if (!MapSheet.isKnown(sheetCode)) throw new Error('Map sheet outside known range: ' + sheetCode); // Shorter tile names for 1:50k if (gridSize === MapSheetTileGridSize) return sheetCode; @@ -548,19 +548,7 @@ export function* iterateMapSheets(bounds: BBox, gridSize: GridSize): Generator Date: Tue, 10 Dec 2024 15:27:43 +1300 Subject: [PATCH 08/12] fix: isDuplicate check --- src/commands/tileindex-validate/tileindex.validate.ts | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/commands/tileindex-validate/tileindex.validate.ts b/src/commands/tileindex-validate/tileindex.validate.ts index f3309a84..90e74e51 100644 --- a/src/commands/tileindex-validate/tileindex.validate.ts +++ b/src/commands/tileindex-validate/tileindex.validate.ts @@ -246,7 +246,7 @@ export const commandTileIndexValidate = command({ return Projection.get(epsg).boundsToGeoJsonFeature(Bounds.fromBbox(loc.bbox), { source: loc.source, tileName: loc.tileNames.join(', '), - isDuplicate: true, // (outputs.get(loc.tileNames)?.length ?? 1) > 1, + isDuplicate: (outputs.get(loc.tileNames.join(', '))?.length ?? 1) > 1, }); }), }); From 641c472fe973e893d0b882ee5a5648f275066ec3 Mon Sep 17 00:00:00 2001 From: Tobias Schmidt <13055656+schmidtnz@users.noreply.github.com> Date: Wed, 11 Dec 2024 12:13:17 +1300 Subject: [PATCH 09/12] fix: getTileName to discard out of bounds sheets --- src/commands/tileindex-validate/tileindex.validate.ts | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/commands/tileindex-validate/tileindex.validate.ts b/src/commands/tileindex-validate/tileindex.validate.ts index 90e74e51..ac90808e 100644 --- a/src/commands/tileindex-validate/tileindex.validate.ts +++ b/src/commands/tileindex-validate/tileindex.validate.ts @@ -502,7 +502,10 @@ export function validateTiffAlignment(tiff: TiffLocation, allowedError = 0.015): export function getTileName(x: number, y: number, gridSize: GridSize): string { const sheetCode = MapSheet.sheetCode(x, y); // TODO: re-enable this check when validation logic - if (!MapSheet.isKnown(sheetCode)) throw new Error('Map sheet outside known range: ' + sheetCode); + if (!MapSheet.isKnown(sheetCode)) { + logger.warn('Map sheet outside known range: ' + sheetCode); + return ''; + } // Shorter tile names for 1:50k if (gridSize === MapSheetTileGridSize) return sheetCode; From 77b1dc664e33e67a8b194033d5b4e681c2cf7f53 Mon Sep 17 00:00:00 2001 From: Tobias Schmidt <13055656+schmidtnz@users.noreply.github.com> Date: Mon, 16 Dec 2024 15:47:44 +1300 Subject: [PATCH 10/12] WIP --- .../__test__/tileindex.validate.test.ts | 59 ++++-- .../tileindex-validate/tileindex.validate.ts | 169 +++++++++++------- 2 files changed, 145 insertions(+), 83 deletions(-) diff --git a/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts b/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts index 4a060dfe..d598b2c1 100644 --- a/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts +++ b/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts @@ -83,23 +83,23 @@ describe('getTileName', () => { describe('tiffLocation', () => { it('get location from tiff', async () => { const TiffAs21 = FakeCogTiff.fromTileName('AS21_1000_0101'); - TiffAs21.images[0].origin[0] = 1492000; - TiffAs21.images[0].origin[1] = 6234000; + // TiffAs21.images[0].origin[0] = 1492000; + // TiffAs21.images[0].origin[1] = 6234000; const TiffAy29 = FakeCogTiff.fromTileName('AY29_1000_0101'); - TiffAy29.images[0].origin[0] = 1684000; - TiffAy29.images[0].origin[1] = 6018000; + // TiffAy29.images[0].origin[0] = 1684000; + // TiffAy29.images[0].origin[1] = 6018000; const location = await extractTiffLocations([TiffAs21, TiffAy29], 1000); - assert.equal(location[0]?.tileNames[0], 'AS21_1000_0101'); - assert.equal(location[1]?.tileNames[0], 'AY29_1000_0101'); + assert.deepEqual(location[0]?.tileNames, ['AS21_1000_0101']); + assert.deepEqual(location[1]?.tileNames, ['AY29_1000_0101']); }); it('should find duplicates', async () => { const TiffAs21 = FakeCogTiff.fromTileName('AS21_1000_0101'); - TiffAs21.images[0].origin[0] = 1492000; - TiffAs21.images[0].origin[1] = 6234000; + // TiffAs21.images[0].origin[0] = 1492000; + // TiffAs21.images[0].origin[1] = 6234000; const TiffAy29 = FakeCogTiff.fromTileName('AY29_1000_0101'); - TiffAy29.images[0].origin[0] = 1684000; - TiffAy29.images[0].origin[1] = 6018000; + // TiffAy29.images[0].origin[0] = 1684000; + // TiffAy29.images[0].origin[1] = 6018000; const location = await extractTiffLocations([TiffAs21, TiffAy29, TiffAs21, TiffAy29], 1000); const duplicates = groupByTileName(location); assert.deepEqual( @@ -118,16 +118,16 @@ describe('tiffLocation', () => { TiffAy29.images[0].origin[0] = 19128043.69337794; TiffAy29.images[0].origin[1] = -4032710.6009459053; const location = await extractTiffLocations([TiffAy29], 1000); - assert.equal(location[0]?.tileNames[0], 'AS21_1000_0101'); + assert.deepEqual(location[0]?.tileNames, ['AS21_1000_0101']); }); it('should fail if one location is not extracted', async () => { const TiffAs21 = FakeCogTiff.fromTileName('AS21_1000_0101'); - TiffAs21.images[0].origin[0] = 1492000; - TiffAs21.images[0].origin[1] = 6234000; + // TiffAs21.images[0].origin[0] = 1492000; + // TiffAs21.images[0].origin[1] = 6234000; const TiffAy29 = FakeCogTiff.fromTileName('AY29_1000_0101'); - TiffAy29.images[0].origin[0] = 1684000; - TiffAy29.images[0].origin[1] = 6018000; + // TiffAy29.images[0].origin[0] = 1684000; + // TiffAy29.images[0].origin[1] = 6018000; TiffAy29.images[0].epsg = 0; // make the projection failing await assert.rejects(extractTiffLocations([TiffAs21, TiffAy29], 1000)); }); @@ -204,7 +204,7 @@ describe('validate', () => { for (const offset of [0.05, -0.05]) { it(`should fail if input tiff origin X is offset by ${offset}m`, async (t) => { - const fakeTiff = FakeCogTiff.fromTileName('AS21_1000_0101'); + const fakeTiff = FakeCogTiff.fromTileName('AU25_1000_0101'); fakeTiff.images[0].origin[0] = fakeTiff.images[0].origin[0] + offset; t.mock.method(TiffLoader, 'load', () => Promise.resolve([fakeTiff])); try { @@ -222,7 +222,7 @@ describe('validate', () => { } }); it(`should fail if input tiff origin Y is offset by ${offset}m`, async (t) => { - const fakeTiff = FakeCogTiff.fromTileName('AS21_1000_0101'); + const fakeTiff = FakeCogTiff.fromTileName('AU25_1000_0101'); fakeTiff.images[0].origin[1] = fakeTiff.images[0].origin[1] + offset; t.mock.method(TiffLoader, 'load', () => Promise.resolve([fakeTiff])); try { @@ -246,7 +246,7 @@ describe('validate', () => { // 720x481 => 720x481 // 721x481 => 721x481 it(`should fail if input tiff width is off by ${offset}m`, async (t) => { - const fakeTiff = FakeCogTiff.fromTileName('AS21_1000_0101'); + const fakeTiff = FakeCogTiff.fromTileName('AU25_1000_0101'); fakeTiff.images[0].size.width = fakeTiff.images[0].size.width + offset; t.mock.method(TiffLoader, 'load', () => Promise.resolve([fakeTiff])); try { @@ -264,7 +264,7 @@ describe('validate', () => { } }); it(`should fail if input tiff height is off by ${offset}m`, async (t) => { - const fakeTiff = FakeCogTiff.fromTileName('AS21_1000_0101'); + const fakeTiff = FakeCogTiff.fromTileName('AU25_1000_0101'); fakeTiff.images[0].size.height = fakeTiff.images[0].size.height + offset; t.mock.method(TiffLoader, 'load', () => Promise.resolve([fakeTiff])); try { @@ -309,3 +309,24 @@ describe('is8BitsTiff', () => { }); }); }); + +describe('TiffFromMisalignedTiff', () => { + it('should properly identify all tiles under a tiff not aligned to our grid', async () => { + const fakeTiffCover4 = FakeCogTiff.fromTileName('CJ09'); + fakeTiffCover4.images[0].origin[0] -= 10; + fakeTiffCover4.images[0].origin[1] += 10; + const fakeTiffCover9 = FakeCogTiff.fromTileName('BA33'); + fakeTiffCover9.images[0].origin[0] -= 10; + fakeTiffCover9.images[0].origin[1] += 10; + fakeTiffCover9.images[0].size.width += 100; + fakeTiffCover9.images[0].size.height += 100; + const fakeTiffCover3 = FakeCogTiff.fromTileName('BL32'); + fakeTiffCover3.images[0].origin[1] -= 10; + fakeTiffCover3.images[0].size.height *= 2; + const locations = await extractTiffLocations([fakeTiffCover4, fakeTiffCover9, fakeTiffCover3], 50000); + + assert.deepEqual(locations[0]?.tileNames, ['CH08', 'CH09', 'CJ08', 'CJ09']); + assert.deepEqual(locations[1]?.tileNames, ['AZ32', 'AZ33', 'AZ34', 'BA32', 'BA33', 'BA34', 'BB32', 'BB33', 'BB34']); + assert.deepEqual(locations[2]?.tileNames, ['BL32', 'BM32', 'BN32']); + }); +}); diff --git a/src/commands/tileindex-validate/tileindex.validate.ts b/src/commands/tileindex-validate/tileindex.validate.ts index ac90808e..d795a101 100644 --- a/src/commands/tileindex-validate/tileindex.validate.ts +++ b/src/commands/tileindex-validate/tileindex.validate.ts @@ -1,5 +1,3 @@ -import assert from 'node:assert'; - import { Bounds, Projection } from '@basemaps/geo'; import { fsa } from '@chunkd/fs'; import { Size, Tiff, TiffTag } from '@cogeotiff/core'; @@ -225,19 +223,23 @@ export const commandTileIndexValidate = command({ } const groupByTileNameStartTime = performance.now(); - const locations = await extractTiffLocations(tiffs, args.scale, args.sourceEpsg); + const tiffLocations = await extractTiffLocations(tiffs, args.scale, args.sourceEpsg); - const outputs = groupByTileName(locations); + const outputTiles = groupByTileName(tiffLocations); logger.info( - { duration: performance.now() - groupByTileNameStartTime, files: locations.length, outputs: outputs.size }, + { + duration: performance.now() - groupByTileNameStartTime, + files: tiffLocations.length, + outputs: outputTiles.size, + }, 'TileIndex: Manifest Assessed for Duplicates', ); if (args.forceOutput || isArgo()) { await fsa.write('/tmp/tile-index-validate/input.geojson', { type: 'FeatureCollection', - features: locations.map((loc) => { + features: tiffLocations.map((loc) => { const epsg = args.sourceEpsg ?? loc.epsg; if (epsg == null) { logger.error({ source: loc.source }, 'TileIndex:Epsg:missing'); @@ -246,7 +248,6 @@ export const commandTileIndexValidate = command({ return Projection.get(epsg).boundsToGeoJsonFeature(Bounds.fromBbox(loc.bbox), { source: loc.source, tileName: loc.tileNames.join(', '), - isDuplicate: (outputs.get(loc.tileNames.join(', '))?.length ?? 1) > 1, }); }), }); @@ -254,11 +255,11 @@ export const commandTileIndexValidate = command({ await fsa.write('/tmp/tile-index-validate/output.geojson', { type: 'FeatureCollection', - features: [...outputs.keys()].map((key) => { + features: [...outputTiles.keys()].map((key) => { const mapTileIndex = MapSheet.getMapTileIndex(key); if (mapTileIndex == null) throw new Error('Failed to extract tile information from: ' + key); return Projection.get(2193).boundsToGeoJsonFeature(Bounds.fromBbox(mapTileIndex.bbox), { - source: outputs.get(key)?.map((l) => l.source), + source: outputTiles.get(key)?.map((l) => l.source), tileName: key, }); }), @@ -267,38 +268,36 @@ export const commandTileIndexValidate = command({ await fsa.write( '/tmp/tile-index-validate/file-list.json', - [...outputs.keys()].map((key) => { - const locs = outputs.get(key); + [...outputTiles.keys()].map((key) => { + const locs = outputTiles.get(key); return { output: key, input: locs?.map((l) => l.source) }; }), ); - logger.info({ path: '/tmp/tile-index-validate/file-list.json', count: outputs.size }, 'Write:FileList'); + logger.info({ path: '/tmp/tile-index-validate/file-list.json', count: outputTiles.size }, 'Write:FileList'); } let retileNeeded = false; - for (const val of outputs.values()) { - if (val.length < 2) continue; - // FIXME + for (const [tileName, tiffs] of outputTiles.entries()) { + if (tiffs.length < 2) continue; if (args.retile) { - const bandType = validateConsistentBands(val); - logger.info( - { tileName: val[0]?.tileNames[0], uris: val.map((v) => v.source), bands: bandType }, - 'TileIndex:Retile', - ); + const bandType = validateConsistentBands(tiffs); + logger.info({ tileName, uris: tiffs.map((v) => v.source), bands: bandType }, 'TileIndex:Retile'); } else { retileNeeded = true; - logger.error({ tileName: val[0]?.tileNames[0], uris: val.map((v) => v.source) }, 'TileIndex:Duplicate'); + logger.error({ tileName, uris: tiffs.map((v) => v.source) }, 'TileIndex:Duplicate'); } } // Validate that all tiffs align to tile grid if (args.validate) { let allValid = true; - for (const tiff of locations) { - const currentValid = validateTiffAlignment(tiff); + for (const tiffLocation of tiffLocations) { + const currentValid = validateTiffAlignment(tiffLocation); allValid = allValid && currentValid; } - if (!allValid) throw new Error(`Tile alignment validation failed`); + if (!allValid) { + throw new Error(`Tile alignment validation failed`); + } } if (retileNeeded) throw new Error(`Duplicate files found, see output.geojson`); @@ -349,7 +348,7 @@ export interface TiffLocation { /** Location to the image */ source: string; /** bbox, [minX, minY, maxX, maxY] */ - bbox: [number, number, number, number]; + bbox: BBox; /** EPSG code of the tiff if found */ epsg?: number | null; /** Output tile name */ @@ -362,6 +361,43 @@ export interface TiffLocation { bands: string[]; } +/** + * Calculate the number grid tiles touched by a TIFF with min being the lowest width/height, and max the largest. + * @param inputMin lowest extent of misaligned tile on one dimension (left or bottom) + * @param inputMax largest extent of misaligned tile on same dimension (right or top) + * @param gridOrigin grid origin (x or y corresponding to min/max dimension) + * @param tileSize height or width of the target grid tile (corresponding to min/max dimension) + */ +const calculateSteps = (inputMin: number, inputMax: number, gridOrigin: number, tileSize: number): number => { + const startStep = Math.floor((inputMin - gridOrigin) / tileSize); + const endStep = Math.floor((inputMax - gridOrigin - 1) / tileSize); + + return endStep - startStep + 1; +}; + +/** + * Reproject the bounding box if the source and target projections are different. + * @param bbox input bounding box + * @param sourceProjection CRS of the input bounding box + * @param targetProjection target CRS + */ +function reprojectIfNeeded( + bbox: BBox, + sourceProjection: Projection, + targetProjection: Projection, +): [number | undefined, number | undefined, number | undefined, number | undefined] { + { + if (targetProjection !== sourceProjection) { + const [ulX, ulY] = targetProjection.fromWgs84(sourceProjection.toWgs84([bbox[0], bbox[3]])); + const [lrX, lrY] = targetProjection.fromWgs84(sourceProjection.toWgs84([bbox[2], bbox[1]])); + + return [ulX, lrY, lrX, ulY]; + } + + return bbox; + } +} + /** * Create a list of `TiffLocation` from a list of TIFFs (`CogTiff`) by extracting their bounding box and generated their tile name from their origin based on a provided `GridSize`. * @@ -386,26 +422,10 @@ export async function extractTiffLocations( return null; } - const centerX = (bbox[0] + bbox[2]) / 2; - const centerY = (bbox[1] + bbox[3]) / 2; - // bbox is not epsg:2193 const targetProjection = Projection.get(2193); const sourceProjection = Projection.get(sourceEpsg); - const [x, y] = targetProjection.fromWgs84(sourceProjection.toWgs84([centerX, centerY])).map(Math.round); - if (x == null || y == null) { - logger.error( - { reason: 'Failed to reproject point', source: tiff.source }, - 'Reprojection:ExtracTiffLocations:Failed', - ); - return null; - } - - // Tilename from center - const tileName = getTileName(x, y, gridSize); - - const [ulX, ulY] = targetProjection.fromWgs84(sourceProjection.toWgs84([bbox[0], bbox[3]])).map(Math.round); - const [lrX, lrY] = targetProjection.fromWgs84(sourceProjection.toWgs84([bbox[2], bbox[1]])).map(Math.round); + const [ulX, lrY, lrX, ulY] = reprojectIfNeeded(bbox, sourceProjection, targetProjection); if (ulX == null || ulY == null || lrX == null || lrY == null) { logger.error( @@ -416,7 +436,6 @@ export async function extractTiffLocations( } const covering = [...iterateMapSheets([ulX, ulY, lrX, lrY], gridSize)] as string[]; - assert.ok(covering.includes(tileName)); // if (shouldValidate) { // // Is the tiff bounding box the same as the map sheet bounding box! @@ -448,13 +467,13 @@ export async function extractTiffLocations( return output; } -export function getSize(extent: [number, number, number, number]): Size { +export function getSize(extent: BBox): Size { return { width: extent[2] - extent[0], height: extent[3] - extent[1] }; } -export function validateTiffAlignment(tiff: TiffLocation, allowedError = 0.015): boolean { +export function validateTiffAlignment(tiff: TiffLocation, allowedErrorMetres = 0.015): boolean { + if (tiff.tileNames.length !== 1) return false; const tileName = tiff.tileNames[0]; - // FIXME if (tileName == null) return false; const mapTileIndex = MapSheet.getMapTileIndex(tileName); if (mapTileIndex == null) { @@ -467,7 +486,7 @@ export function validateTiffAlignment(tiff: TiffLocation, allowedError = 0.015): // Top Left const errX = Math.abs(tiff.bbox[0] - mapTileIndex.bbox[0]); const errY = Math.abs(tiff.bbox[3] - mapTileIndex.bbox[3]); - if (errX > allowedError || errY > allowedError) { + if (errX > allowedErrorMetres || errY > allowedErrorMetres) { logger.error( { reason: `The origin is invalid x:${tiff.bbox[0]}, y:${tiff.bbox[3]}`, source: tiff.source }, 'TileInvalid:Validation:Failed', @@ -475,7 +494,6 @@ export function validateTiffAlignment(tiff: TiffLocation, allowedError = 0.015): return false; } - // TODO do we validate bottom right const tiffSize = getSize(tiff.bbox); if (tiffSize.width !== mapTileIndex.width) { logger.error( @@ -501,10 +519,8 @@ export function validateTiffAlignment(tiff: TiffLocation, allowedError = 0.015): export function getTileName(x: number, y: number, gridSize: GridSize): string { const sheetCode = MapSheet.sheetCode(x, y); - // TODO: re-enable this check when validation logic if (!MapSheet.isKnown(sheetCode)) { - logger.warn('Map sheet outside known range: ' + sheetCode); - return ''; + logger.info(`Map sheet (${sheetCode}) at coordinates (${x}, ${y}) is outside the known range.`); } // Shorter tile names for 1:50k @@ -516,8 +532,8 @@ export function getTileName(x: number, y: number, gridSize: GridSize): string { const nbDigits = gridSize === 500 ? 3 : 2; - const offsetX = Math.round(Math.floor((x - MapSheet.origin.x) / MapSheet.width)); - const offsetY = Math.round(Math.floor((MapSheet.origin.y - y) / MapSheet.height)); + const offsetX = Math.floor((x - MapSheet.origin.x) / MapSheet.width); + const offsetY = Math.floor((MapSheet.origin.y - y) / MapSheet.height); const maxY = MapSheet.origin.y - offsetY * MapSheet.height; const minX = MapSheet.origin.x + offsetX * MapSheet.width; const tileX = Math.round(Math.floor((x - minX) / tileWidth + 1)); @@ -545,19 +561,44 @@ export async function validate8BitsTiff(tiff: Tiff): Promise { } } -export function* iterateMapSheets(bounds: BBox, gridSize: GridSize): Generator { - const minX = Math.min(bounds[0], bounds[2]); - const maxX = Math.max(bounds[0], bounds[2]); - const minY = Math.min(bounds[1], bounds[3]); - const maxY = Math.max(bounds[1], bounds[3]); +export function* iterateMapSheets(bbox: BBox, gridSize: GridSize): Generator { + const bounds = Bounds.fromBbox(bbox); - const tilesPerMapSheet = Math.floor(MapSheet.gridSizeMax / gridSize); + const minX = Math.min(bbox[0], bbox[2]); + const maxX = Math.max(bbox[0], bbox[2]); + const minY = Math.min(bbox[1], bbox[3]); + const maxY = Math.max(bbox[1], bbox[3]); - const tileWidth = Math.floor(MapSheet.width / tilesPerMapSheet); - const tileHeight = Math.floor(MapSheet.height / tilesPerMapSheet); - for (let x = minX; x < maxX; x += tileWidth) { - for (let y = maxY; y > minY; y -= tileHeight) { - yield getTileName(x, y, gridSize); + const tilesPerMapSheetSide = Math.floor(MapSheet.gridSizeMax / gridSize); + + const tileWidth = Math.floor(MapSheet.width / tilesPerMapSheetSide); + const tileHeight = Math.floor(MapSheet.height / tilesPerMapSheetSide); + + const stepsX = calculateSteps(minX, maxX, MapSheet.origin.x, tileWidth); + const stepsY = calculateSteps(-maxY, -minY, -MapSheet.origin.y, tileHeight); // Inverted Y axis + + for (let stepY = 0; stepY < stepsY; stepY += 1) { + const y = maxY - stepY * tileHeight; + for (let stepX = 0; stepX < stepsX; stepX += 1) { + const x = minX + stepX * tileWidth; + const tileName = getTileName(x, y, gridSize); + const tile = MapSheet.getMapTileIndex(tileName); + if (tile == null) { + logger.error({ tileName, x, y, gridSize }, `No MapTileIndex for tile ${tileName} at xy ${x}, ${y}`); + continue; + } + const intersection = bounds.intersection(Bounds.fromBbox(tile.bbox)); + if (intersection === null) { + logger.warn({ tileName, tile, x, y, gridSize }, `Source TIFF at xy ${x}, ${y} does not intersect tile. Skipping ${tileName}`); + } else if (intersection?.width <= 1 || intersection?.height <= 1) { + // TODO: Check if we can make this relative to GSD and use <1 pixel instead of meters. + logger.warn( + { tileName, tile, x, y, gridSize, intersection }, + `Minor intersection w${intersection.width} by h${intersection.height} at xy ${x}, ${y}. Skipping tile ${tileName}`, + ); + } else { + yield tileName; + } } } } From 94ed7b722aa37dc764220463796771383924c68b Mon Sep 17 00:00:00 2001 From: Tobias Schmidt <13055656+schmidtnz@users.noreply.github.com> Date: Thu, 19 Dec 2024 12:06:54 +1300 Subject: [PATCH 11/12] feat: retile misaligned input TIFFs --- .../__test__/tileindex.validate.test.ts | 60 +++++--- .../tileindex-validate/tileindex.validate.ts | 132 +++++++++--------- 2 files changed, 108 insertions(+), 84 deletions(-) diff --git a/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts b/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts index d598b2c1..52812dc8 100644 --- a/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts +++ b/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts @@ -1,8 +1,10 @@ import assert from 'node:assert'; import { before, beforeEach, describe, it } from 'node:test'; +import { Projection } from '@basemaps/geo'; import { fsa } from '@chunkd/fs'; import { FsMemory } from '@chunkd/source-memory'; +import { BBox } from '@linzjs/geojson'; import { FeatureCollection } from 'geojson'; import { MapSheetData } from '../../../utils/__test__/mapsheet.data.js'; @@ -14,6 +16,7 @@ import { getTileName, GridSizeFromString, groupByTileName, + reprojectIfNeeded, TiffLoader, validate8BitsTiff, } from '../tileindex.validate.js'; @@ -83,11 +86,11 @@ describe('getTileName', () => { describe('tiffLocation', () => { it('get location from tiff', async () => { const TiffAs21 = FakeCogTiff.fromTileName('AS21_1000_0101'); - // TiffAs21.images[0].origin[0] = 1492000; - // TiffAs21.images[0].origin[1] = 6234000; + TiffAs21.images[0].origin[0] = 1492000; + TiffAs21.images[0].origin[1] = 6234000; const TiffAy29 = FakeCogTiff.fromTileName('AY29_1000_0101'); - // TiffAy29.images[0].origin[0] = 1684000; - // TiffAy29.images[0].origin[1] = 6018000; + TiffAy29.images[0].origin[0] = 1684000; + TiffAy29.images[0].origin[1] = 6018000; const location = await extractTiffLocations([TiffAs21, TiffAy29], 1000); assert.deepEqual(location[0]?.tileNames, ['AS21_1000_0101']); assert.deepEqual(location[1]?.tileNames, ['AY29_1000_0101']); @@ -95,11 +98,11 @@ describe('tiffLocation', () => { it('should find duplicates', async () => { const TiffAs21 = FakeCogTiff.fromTileName('AS21_1000_0101'); - // TiffAs21.images[0].origin[0] = 1492000; - // TiffAs21.images[0].origin[1] = 6234000; + TiffAs21.images[0].origin[0] = 1492000; + TiffAs21.images[0].origin[1] = 6234000; const TiffAy29 = FakeCogTiff.fromTileName('AY29_1000_0101'); - // TiffAy29.images[0].origin[0] = 1684000; - // TiffAy29.images[0].origin[1] = 6018000; + TiffAy29.images[0].origin[0] = 1684000; + TiffAy29.images[0].origin[1] = 6018000; const location = await extractTiffLocations([TiffAs21, TiffAy29, TiffAs21, TiffAy29], 1000); const duplicates = groupByTileName(location); assert.deepEqual( @@ -123,11 +126,11 @@ describe('tiffLocation', () => { it('should fail if one location is not extracted', async () => { const TiffAs21 = FakeCogTiff.fromTileName('AS21_1000_0101'); - // TiffAs21.images[0].origin[0] = 1492000; - // TiffAs21.images[0].origin[1] = 6234000; + TiffAs21.images[0].origin[0] = 1492000; + TiffAs21.images[0].origin[1] = 6234000; const TiffAy29 = FakeCogTiff.fromTileName('AY29_1000_0101'); - // TiffAy29.images[0].origin[0] = 1684000; - // TiffAy29.images[0].origin[1] = 6018000; + TiffAy29.images[0].origin[0] = 1684000; + TiffAy29.images[0].origin[1] = 6018000; TiffAy29.images[0].epsg = 0; // make the projection failing await assert.rejects(extractTiffLocations([TiffAs21, TiffAy29], 1000)); }); @@ -204,7 +207,7 @@ describe('validate', () => { for (const offset of [0.05, -0.05]) { it(`should fail if input tiff origin X is offset by ${offset}m`, async (t) => { - const fakeTiff = FakeCogTiff.fromTileName('AU25_1000_0101'); + const fakeTiff = FakeCogTiff.fromTileName('AS21_1000_0101'); fakeTiff.images[0].origin[0] = fakeTiff.images[0].origin[0] + offset; t.mock.method(TiffLoader, 'load', () => Promise.resolve([fakeTiff])); try { @@ -222,7 +225,7 @@ describe('validate', () => { } }); it(`should fail if input tiff origin Y is offset by ${offset}m`, async (t) => { - const fakeTiff = FakeCogTiff.fromTileName('AU25_1000_0101'); + const fakeTiff = FakeCogTiff.fromTileName('AS21_1000_0101'); fakeTiff.images[0].origin[1] = fakeTiff.images[0].origin[1] + offset; t.mock.method(TiffLoader, 'load', () => Promise.resolve([fakeTiff])); try { @@ -246,7 +249,7 @@ describe('validate', () => { // 720x481 => 720x481 // 721x481 => 721x481 it(`should fail if input tiff width is off by ${offset}m`, async (t) => { - const fakeTiff = FakeCogTiff.fromTileName('AU25_1000_0101'); + const fakeTiff = FakeCogTiff.fromTileName('AS21_1000_0101'); fakeTiff.images[0].size.width = fakeTiff.images[0].size.width + offset; t.mock.method(TiffLoader, 'load', () => Promise.resolve([fakeTiff])); try { @@ -264,7 +267,7 @@ describe('validate', () => { } }); it(`should fail if input tiff height is off by ${offset}m`, async (t) => { - const fakeTiff = FakeCogTiff.fromTileName('AU25_1000_0101'); + const fakeTiff = FakeCogTiff.fromTileName('AS21_1000_0101'); fakeTiff.images[0].size.height = fakeTiff.images[0].size.height + offset; t.mock.method(TiffLoader, 'load', () => Promise.resolve([fakeTiff])); try { @@ -330,3 +333,28 @@ describe('TiffFromMisalignedTiff', () => { assert.deepEqual(locations[2]?.tileNames, ['BL32', 'BM32', 'BN32']); }); }); + +describe('reprojectIfNeeded', () => { + it('should reproject the bounding box if projections are different', () => { + const sourceProjection = Projection.get(4326); // WGS84 + const targetProjection = Projection.get(2193); // New Zealand Transverse Mercator 2000 + const bbox: BBox = [172, -41, 174, -38]; // Example bounding box in WGS84 + + const reprojectedBbox = reprojectIfNeeded(bbox, sourceProjection, targetProjection) as BBox; + assert.notDeepEqual(bbox.map(Math.round), reprojectedBbox.map(Math.round)); // expect output to be quite different + assert.ok(reprojectedBbox); + + const roundtripBbox = reprojectIfNeeded(reprojectedBbox, targetProjection, sourceProjection); + assert.deepEqual(bbox.map(Math.round), roundtripBbox?.map(Math.round)); // expect output to be very similar (floating point / rounding errors) + }); + + it('should return the same bounding box if projections are the same', () => { + const sourceProjection = Projection.get(2193); // New Zealand Transverse Mercator 2000 + const targetProjection = Projection.get(2193); // New Zealand Transverse Mercator 2000 + const bbox: BBox = [1, 2, 3, 4]; // Example bounding box + + const reprojectedBbox = reprojectIfNeeded(bbox, sourceProjection, targetProjection); + + assert.deepEqual(reprojectedBbox, bbox); + }); +}); diff --git a/src/commands/tileindex-validate/tileindex.validate.ts b/src/commands/tileindex-validate/tileindex.validate.ts index d795a101..26af9623 100644 --- a/src/commands/tileindex-validate/tileindex.validate.ts +++ b/src/commands/tileindex-validate/tileindex.validate.ts @@ -278,7 +278,8 @@ export const commandTileIndexValidate = command({ let retileNeeded = false; for (const [tileName, tiffs] of outputTiles.entries()) { - if (tiffs.length < 2) continue; + if (tiffs.length === 0) throw new Error(`Output tile with no source tiff: ${tileName}`); + if (tiffs.length === 1) continue; if (args.retile) { const bandType = validateConsistentBands(tiffs); logger.info({ tileName, uris: tiffs.map((v) => v.source), bands: bandType }, 'TileIndex:Retile'); @@ -361,40 +362,18 @@ export interface TiffLocation { bands: string[]; } -/** - * Calculate the number grid tiles touched by a TIFF with min being the lowest width/height, and max the largest. - * @param inputMin lowest extent of misaligned tile on one dimension (left or bottom) - * @param inputMax largest extent of misaligned tile on same dimension (right or top) - * @param gridOrigin grid origin (x or y corresponding to min/max dimension) - * @param tileSize height or width of the target grid tile (corresponding to min/max dimension) - */ -const calculateSteps = (inputMin: number, inputMax: number, gridOrigin: number, tileSize: number): number => { - const startStep = Math.floor((inputMin - gridOrigin) / tileSize); - const endStep = Math.floor((inputMax - gridOrigin - 1) / tileSize); - - return endStep - startStep + 1; -}; - /** * Reproject the bounding box if the source and target projections are different. * @param bbox input bounding box * @param sourceProjection CRS of the input bounding box * @param targetProjection target CRS */ -function reprojectIfNeeded( - bbox: BBox, - sourceProjection: Projection, - targetProjection: Projection, -): [number | undefined, number | undefined, number | undefined, number | undefined] { +export function reprojectIfNeeded(bbox: BBox, sourceProjection: Projection, targetProjection: Projection): BBox | null { { - if (targetProjection !== sourceProjection) { - const [ulX, ulY] = targetProjection.fromWgs84(sourceProjection.toWgs84([bbox[0], bbox[3]])); - const [lrX, lrY] = targetProjection.fromWgs84(sourceProjection.toWgs84([bbox[2], bbox[1]])); - - return [ulX, lrY, lrX, ulY]; - } - - return bbox; + if (targetProjection === sourceProjection) return bbox; + const [ulX, ulY] = targetProjection.fromWgs84(sourceProjection.toWgs84([bbox[0], bbox[3]])) as [number, number]; // Typescript inferred this as number[] but should be [number, number] | [number, number, number] + const [lrX, lrY] = targetProjection.fromWgs84(sourceProjection.toWgs84([bbox[2], bbox[1]])) as [number, number]; // Typescript inferred this as number[] but should be [number, number] | [number, number, number] + return [Math.min(ulX, lrX), Math.min(lrY, ulY), Math.max(ulX, lrX), Math.max(lrY, ulY)]; } } @@ -414,7 +393,7 @@ export async function extractTiffLocations( const result = await Promise.all( tiffs.map(async (tiff): Promise => { try { - const bbox = await findBoundingBox(tiff); + const sourceBbox = await findBoundingBox(tiff); const sourceEpsg = forceSourceEpsg ?? tiff.images[0]?.epsg; if (sourceEpsg == null) { @@ -425,9 +404,9 @@ export async function extractTiffLocations( const targetProjection = Projection.get(2193); const sourceProjection = Projection.get(sourceEpsg); - const [ulX, lrY, lrX, ulY] = reprojectIfNeeded(bbox, sourceProjection, targetProjection); + const targetBbox = reprojectIfNeeded(sourceBbox, sourceProjection, targetProjection); - if (ulX == null || ulY == null || lrX == null || lrY == null) { + if (targetBbox === null) { logger.error( { reason: 'Failed to reproject point', source: tiff.source }, 'Reprojection:ExtracTiffLocations:Failed', @@ -435,7 +414,7 @@ export async function extractTiffLocations( return null; } - const covering = [...iterateMapSheets([ulX, ulY, lrX, lrY], gridSize)] as string[]; + const covering = getCovering(targetBbox, gridSize); // if (shouldValidate) { // // Is the tiff bounding box the same as the map sheet bounding box! @@ -443,7 +422,7 @@ export async function extractTiffLocations( // // assert bbox == MapSheet.getMapTileIndex(tileName).bbox // } return { - bbox, + bbox: targetBbox, source: tiff.source.url.href, tileNames: covering, epsg: tiff.images[0]?.epsg, @@ -520,7 +499,10 @@ export function validateTiffAlignment(tiff: TiffLocation, allowedErrorMetres = 0 export function getTileName(x: number, y: number, gridSize: GridSize): string { const sheetCode = MapSheet.sheetCode(x, y); if (!MapSheet.isKnown(sheetCode)) { - logger.info(`Map sheet (${sheetCode}) at coordinates (${x}, ${y}) is outside the known range.`); + logger.info( + { sheetCode, x, y, gridSize }, + `Map sheet (${sheetCode}) at coordinates (${x}, ${y}) is outside the known range.`, + ); } // Shorter tile names for 1:50k @@ -561,44 +543,58 @@ export async function validate8BitsTiff(tiff: Tiff): Promise { } } -export function* iterateMapSheets(bbox: BBox, gridSize: GridSize): Generator { - const bounds = Bounds.fromBbox(bbox); - - const minX = Math.min(bbox[0], bbox[2]); - const maxX = Math.max(bbox[0], bbox[2]); - const minY = Math.min(bbox[1], bbox[3]); - const maxY = Math.max(bbox[1], bbox[3]); - +/** + * Get the list of map sheets / tiles that intersect with the given bounding box. + * + * @param bbox Bounding box of the area of interest (in EPSG:2193) to get the map sheets for (e.g. TIFF area). + * @param gridSize Grid size of the map sheets / tiles to get. + * @param minIntersectionMeters Minimum intersection area in meters (width or height) to include the map sheet. + */ +function getCovering(bbox: BBox, gridSize: GridSize, minIntersectionMeters = 0.15): string[] { + const SurroundingTiles = [ + { x: 1, y: 0 }, + { x: 0, y: -1 }, // inverted Y axis + // Future versions may want to explore tiles in all directions + // { x: 0, y: 1 }, + // { x: -1, y: 0 }, + ]; + + const targetBounds = Bounds.fromBbox(bbox); + + const output: string[] = []; const tilesPerMapSheetSide = Math.floor(MapSheet.gridSizeMax / gridSize); const tileWidth = Math.floor(MapSheet.width / tilesPerMapSheetSide); const tileHeight = Math.floor(MapSheet.height / tilesPerMapSheetSide); - const stepsX = calculateSteps(minX, maxX, MapSheet.origin.x, tileWidth); - const stepsY = calculateSteps(-maxY, -minY, -MapSheet.origin.y, tileHeight); // Inverted Y axis - - for (let stepY = 0; stepY < stepsY; stepY += 1) { - const y = maxY - stepY * tileHeight; - for (let stepX = 0; stepX < stepsX; stepX += 1) { - const x = minX + stepX * tileWidth; - const tileName = getTileName(x, y, gridSize); - const tile = MapSheet.getMapTileIndex(tileName); - if (tile == null) { - logger.error({ tileName, x, y, gridSize }, `No MapTileIndex for tile ${tileName} at xy ${x}, ${y}`); - continue; - } - const intersection = bounds.intersection(Bounds.fromBbox(tile.bbox)); - if (intersection === null) { - logger.warn({ tileName, tile, x, y, gridSize }, `Source TIFF at xy ${x}, ${y} does not intersect tile. Skipping ${tileName}`); - } else if (intersection?.width <= 1 || intersection?.height <= 1) { - // TODO: Check if we can make this relative to GSD and use <1 pixel instead of meters. - logger.warn( - { tileName, tile, x, y, gridSize, intersection }, - `Minor intersection w${intersection.width} by h${intersection.height} at xy ${x}, ${y}. Skipping tile ${tileName}`, - ); - } else { - yield tileName; - } - } + const seen = new Set(); + const todo: Bounds[] = []; + + const sheetName = getTileName(bbox[0], bbox[3], gridSize); + const sheetInfo = MapSheet.getMapTileIndex(sheetName); + if (sheetInfo == null) throw new Error('Unable to extract sheet information for point: ' + bbox[0]); + + todo.push(Bounds.fromBbox(sheetInfo.bbox)); + while (todo.length > 0) { + const nextBounds = todo.shift(); + if (nextBounds == null) continue; + + const nextX = nextBounds.x; + const nextY = nextBounds.bottom; // inverted Y axis + const bboxId = `${nextX}:${nextY}:${nextBounds.width}:${nextBounds.height}`; + // Only process each sheet once + if (seen.has(bboxId)) continue; + seen.add(bboxId); + + const intersection = targetBounds.intersection(nextBounds); + if (intersection == null) continue; // no intersection, target mapshet is outside bounds of source + // Check all the surrounding tiles + for (const pt of SurroundingTiles) todo.push(nextBounds.add({ x: pt.x * tileWidth, y: pt.y * tileHeight })); // intersection not null, so add all neighbours + // Add to output only if the intersection is above the minimum coverage + if (intersection.width < minIntersectionMeters || intersection.height < minIntersectionMeters) continue; + output.push(getTileName(nextX, nextY, gridSize)); } + + output.sort(); + return output; } From afad58da70baf2d87dc0dacf7b86306a951db013 Mon Sep 17 00:00:00 2001 From: Tobias Schmidt <13055656+schmidtnz@users.noreply.github.com> Date: Fri, 20 Dec 2024 15:58:30 +1300 Subject: [PATCH 12/12] feat: retile misaligned input TIFFs --- .../tileindex-validate/tileindex.validate.ts | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/commands/tileindex-validate/tileindex.validate.ts b/src/commands/tileindex-validate/tileindex.validate.ts index 26af9623..72ed83f1 100644 --- a/src/commands/tileindex-validate/tileindex.validate.ts +++ b/src/commands/tileindex-validate/tileindex.validate.ts @@ -371,8 +371,10 @@ export interface TiffLocation { export function reprojectIfNeeded(bbox: BBox, sourceProjection: Projection, targetProjection: Projection): BBox | null { { if (targetProjection === sourceProjection) return bbox; - const [ulX, ulY] = targetProjection.fromWgs84(sourceProjection.toWgs84([bbox[0], bbox[3]])) as [number, number]; // Typescript inferred this as number[] but should be [number, number] | [number, number, number] - const [lrX, lrY] = targetProjection.fromWgs84(sourceProjection.toWgs84([bbox[2], bbox[1]])) as [number, number]; // Typescript inferred this as number[] but should be [number, number] | [number, number, number] + // fromWgs84 and toWgs84 functions are typed as number[] but could be refined to [number, number] | [number, number, number]. + // With 2 input args, they will return [number, number]. + const [ulX, ulY] = targetProjection.fromWgs84(sourceProjection.toWgs84([bbox[0], bbox[3]])) as [number, number]; + const [lrX, lrY] = targetProjection.fromWgs84(sourceProjection.toWgs84([bbox[2], bbox[1]])) as [number, number]; return [Math.min(ulX, lrX), Math.min(lrY, ulY), Math.max(ulX, lrX), Math.max(lrY, ulY)]; } } @@ -553,10 +555,9 @@ export async function validate8BitsTiff(tiff: Tiff): Promise { function getCovering(bbox: BBox, gridSize: GridSize, minIntersectionMeters = 0.15): string[] { const SurroundingTiles = [ { x: 1, y: 0 }, - { x: 0, y: -1 }, // inverted Y axis - // Future versions may want to explore tiles in all directions - // { x: 0, y: 1 }, - // { x: -1, y: 0 }, + { x: 0, y: -1 }, + { x: -1, y: 0 }, + { x: 0, y: 1 }, ]; const targetBounds = Bounds.fromBbox(bbox); @@ -595,6 +596,5 @@ function getCovering(bbox: BBox, gridSize: GridSize, minIntersectionMeters = 0.1 output.push(getTileName(nextX, nextY, gridSize)); } - output.sort(); - return output; + return output.sort(); }