diff --git a/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts b/src/commands/tileindex-validate/__test__/tileindex.validate.test.ts index f182c88d..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'; @@ -89,8 +92,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.deepEqual(location[0]?.tileNames, ['AS21_1000_0101']); + assert.deepEqual(location[1]?.tileNames, ['AY29_1000_0101']); }); it('should find duplicates', async () => { @@ -118,7 +121,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.deepEqual(location[0]?.tileNames, ['AS21_1000_0101']); }); it('should fail if one location is not extracted', async () => { @@ -309,3 +312,49 @@ 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']); + }); +}); + +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 f97c5e5b..72ed83f1 100644 --- a/src/commands/tileindex-validate/tileindex.validate.ts +++ b/src/commands/tileindex-validate/tileindex.validate.ts @@ -1,6 +1,7 @@ import { Bounds, 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'; @@ -9,7 +10,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'; @@ -222,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'); @@ -242,8 +247,7 @@ 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(', '), }); }), }); @@ -251,14 +255,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: [...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: locs.map((l) => l.source), - tileName: firstLoc.tileName, + source: outputTiles.get(key)?.map((l) => l.source), + tileName: key, }); }), }); @@ -266,36 +268,37 @@ 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) }; + [...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; + for (const [tileName, tiffs] of outputTiles.entries()) { + 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(val); - logger.info( - { tileName: val[0]?.tileName, 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]?.tileName, 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`); @@ -333,9 +336,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; } @@ -344,11 +349,11 @@ 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 */ - tileName: string; + tileNames: string[]; /** * List of bands inside the tiff in the format `uint8` `uint16` * @@ -357,6 +362,23 @@ export interface TiffLocation { bands: string[]; } +/** + * 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 + */ +export function reprojectIfNeeded(bbox: BBox, sourceProjection: Projection, targetProjection: Projection): BBox | null { + { + if (targetProjection === sourceProjection) return bbox; + // 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)]; + } +} + /** * 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`. * @@ -373,7 +395,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) { @@ -381,14 +403,12 @@ 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])); - if (x == null || y == null) { + const targetBbox = reprojectIfNeeded(sourceBbox, sourceProjection, targetProjection); + + if (targetBbox === null) { logger.error( { reason: 'Failed to reproject point', source: tiff.source }, 'Reprojection:ExtracTiffLocations:Failed', @@ -396,8 +416,7 @@ export async function extractTiffLocations( return null; } - // Tilename from center - const tileName = getTileName(x, y, gridSize); + const covering = getCovering(targetBbox, gridSize); // if (shouldValidate) { // // Is the tiff bounding box the same as the map sheet bounding box! @@ -405,9 +424,9 @@ export async function extractTiffLocations( // // assert bbox == MapSheet.getMapTileIndex(tileName).bbox // } return { - bbox, + bbox: targetBbox, source: tiff.source.url.href, - tileName, + tileNames: covering, epsg: tiff.images[0]?.epsg, bands: await extractBandInformation(tiff), }; @@ -429,15 +448,18 @@ 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 { - const mapTileIndex = MapSheet.getMapTileIndex(tiff.tileName); +export function validateTiffAlignment(tiff: TiffLocation, allowedErrorMetres = 0.015): boolean { + if (tiff.tileNames.length !== 1) return false; + const tileName = tiff.tileNames[0]; + 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; @@ -445,7 +467,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', @@ -453,7 +475,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( @@ -478,14 +499,13 @@ 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')}`; - // TODO: re-enable this check when validation logic - // if (!MapSheet.isKnown(sheetCode)) throw new Error('Map sheet outside known range: ' + sheetCode); + const sheetCode = MapSheet.sheetCode(x, y); + if (!MapSheet.isKnown(sheetCode)) { + logger.info( + { sheetCode, x, y, gridSize }, + `Map sheet (${sheetCode}) at coordinates (${x}, ${y}) is outside the known range.`, + ); + } // Shorter tile names for 1:50k if (gridSize === MapSheetTileGridSize) return sheetCode; @@ -496,6 +516,8 @@ export function getTileName(x: number, y: number, gridSize: GridSize): string { const nbDigits = gridSize === 500 ? 3 : 2; + 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)); @@ -522,3 +544,57 @@ export async function validate8BitsTiff(tiff: Tiff): Promise { throw new Error(`${tiff.source.url.href} is not a 8 bits TIFF`); } } + +/** + * 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 }, + { x: -1, y: 0 }, + { x: 0, y: 1 }, + ]; + + 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 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)); + } + + return output.sort(); +}