From f799006ccf1a45ec8aebfe132603a17c031e4825 Mon Sep 17 00:00:00 2001 From: Geoff Jacobsen Date: Fri, 3 Jul 2020 16:08:55 +1200 Subject: [PATCH] fix(cli): refactor projection logic to allow chathams to be built (#854) Source projection does not need to have a TileMatrixSet. This commit removes that constraint. --- .../cli/cogify/__test__/action.batch.test.ts | 2 +- packages/cli/src/cli/cogify/action.batch.ts | 4 +- packages/cli/src/cog/__test__/cog.test.ts | 2 +- packages/cli/src/cog/__test__/cog.vrt.test.ts | 20 +-- packages/cli/src/cog/__test__/cutline.test.ts | 22 ++-- packages/cli/src/cog/builder.ts | 12 +- packages/cli/src/cog/cog.ts | 12 +- packages/cli/src/cog/cog.vrt.ts | 2 +- packages/cli/src/cog/cutline.ts | 40 +++--- packages/cli/src/cog/job.ts | 11 +- .../aws/__test__/tile.metadata.table.test.ts | 2 +- packages/shared/src/index.ts | 2 +- .../src/proj/__test__/projection.test.ts | 68 ++++++++++ .../projection.tile.matrix.set.test.ts | 107 ++++------------ .../src/{tms => proj}/__test__/test.util.ts | 4 +- packages/shared/src/{tms => proj}/citm2000.ts | 0 packages/shared/src/{tms => proj}/nztm2000.ts | 0 .../projection.tile.matrix.set.ts | 79 +----------- packages/shared/src/proj/projection.ts | 116 ++++++++++++++++++ 19 files changed, 286 insertions(+), 219 deletions(-) create mode 100644 packages/shared/src/proj/__test__/projection.test.ts rename packages/shared/src/{tms => proj}/__test__/projection.tile.matrix.set.test.ts (71%) rename packages/shared/src/{tms => proj}/__test__/test.util.ts (85%) rename packages/shared/src/{tms => proj}/citm2000.ts (100%) rename packages/shared/src/{tms => proj}/nztm2000.ts (100%) rename packages/shared/src/{tms => proj}/projection.tile.matrix.set.ts (57%) create mode 100644 packages/shared/src/proj/projection.ts diff --git a/packages/cli/src/cli/cogify/__test__/action.batch.test.ts b/packages/cli/src/cli/cogify/__test__/action.batch.test.ts index 15af5ffc5..e75de199f 100644 --- a/packages/cli/src/cli/cogify/__test__/action.batch.test.ts +++ b/packages/cli/src/cli/cogify/__test__/action.batch.test.ts @@ -1,7 +1,7 @@ import o from 'ospec'; import { CogJob } from '../../../cog/types'; import { extractResolutionFromName, extractYearFromName, createImageryRecordFromJob } from '../action.batch'; -import { qkToNamedBounds } from '@basemaps/shared/build/tms/__test__/test.util'; +import { qkToNamedBounds } from '@basemaps/shared/build/proj/__test__/test.util'; import { Bounds } from '@basemaps/geo'; import { round } from '@basemaps/test/build/rounding'; diff --git a/packages/cli/src/cli/cogify/action.batch.ts b/packages/cli/src/cli/cogify/action.batch.ts index 4efd2e376..e25dc3d57 100644 --- a/packages/cli/src/cli/cogify/action.batch.ts +++ b/packages/cli/src/cli/cogify/action.batch.ts @@ -103,9 +103,9 @@ export class ActionBatchJob extends CommandLineAction { isCommit: boolean, ): Promise<{ jobName: string; jobId: string; memory: number }> { const jobName = ActionBatchJob.id(job, name); - const targetProj = ProjectionTileMatrixSet.get(job.projection); + const targetPtms = ProjectionTileMatrixSet.get(job.projection); const tile = TileMatrixSet.nameToTile(name); - const alignmentLevels = targetProj.findAlignmentLevels(tile, job.source.pixelScale); + const alignmentLevels = targetPtms.findAlignmentLevels(tile, job.source.pixelScale); // Give 25% more memory to larger jobs const resDiff = 1 + Math.max(alignmentLevels - MagicAlignmentLevel, 0) * 0.25; const memory = 3900 * resDiff; diff --git a/packages/cli/src/cog/__test__/cog.test.ts b/packages/cli/src/cog/__test__/cog.test.ts index ee8b71bad..5f33ae20b 100644 --- a/packages/cli/src/cog/__test__/cog.test.ts +++ b/packages/cli/src/cog/__test__/cog.test.ts @@ -6,7 +6,7 @@ import { buildCogForName } from '../cog'; import { SourceTiffTestHelper } from './source.tiff.testhelper'; import { TilingScheme } from '../../gdal/gdal.config'; import { round } from '@basemaps/test/build/rounding'; -import { qkToName } from '@basemaps/shared/build/tms/__test__/test.util'; +import { qkToName } from '@basemaps/shared/build/proj/__test__/test.util'; LogConfig.disable(); diff --git a/packages/cli/src/cog/__test__/cog.vrt.test.ts b/packages/cli/src/cog/__test__/cog.vrt.test.ts index 9a1003ccb..d87343a66 100644 --- a/packages/cli/src/cog/__test__/cog.vrt.test.ts +++ b/packages/cli/src/cog/__test__/cog.vrt.test.ts @@ -1,6 +1,6 @@ import { EpsgCode } from '@basemaps/geo'; import { FileOperatorSimple, LogConfig, ProjectionTileMatrixSet } from '@basemaps/shared'; -import { qkToName } from '@basemaps/shared/build/tms/__test__/test.util'; +import { qkToName } from '@basemaps/shared/build/proj/__test__/test.util'; import { round } from '@basemaps/test/build/rounding'; import o from 'ospec'; import { GdalCogBuilder } from '../../gdal/gdal'; @@ -18,7 +18,7 @@ o.spec('cog.vrt', () => { const testDir = `${__dirname}/../../../__test.assets__`; - const googleProj = ProjectionTileMatrixSet.get(EpsgCode.Google); + const googlePtms = ProjectionTileMatrixSet.get(EpsgCode.Google); const sourceBounds = SourceTiffTestHelper.tiffNztmBounds(testDir); const [tif1, tif2] = sourceBounds; @@ -56,8 +56,8 @@ o.spec('cog.vrt', () => { }); o('1 crosses, 1 outside', async () => { - const cutline = new Cutline(googleProj, await Cutline.loadCutline(testDir + '/kapiti.geojson')); - const cl2 = new Cutline(googleProj, await Cutline.loadCutline(testDir + '/mana.geojson')); + const cutline = new Cutline(googlePtms, await Cutline.loadCutline(testDir + '/kapiti.geojson')); + const cl2 = new Cutline(googlePtms, await Cutline.loadCutline(testDir + '/mana.geojson')); cutline.clipPoly.push(...cl2.clipPoly); job.source.resZoom = 17; @@ -70,7 +70,7 @@ o.spec('cog.vrt', () => { }); o('not within tile', async () => { - const cutline = new Cutline(googleProj, await Cutline.loadCutline(testDir + '/kapiti.geojson')); + const cutline = new Cutline(googlePtms, await Cutline.loadCutline(testDir + '/kapiti.geojson')); const vrt = await CogVrt.buildVrt(tmpFolder, job, cutline, qkToName('3131110001'), logger); @@ -80,7 +80,7 @@ o.spec('cog.vrt', () => { }); o('no cutline same projection', async () => { - const vrt = await CogVrt.buildVrt(tmpFolder, job, new Cutline(googleProj), qkToName('31'), logger); + const vrt = await CogVrt.buildVrt(tmpFolder, job, new Cutline(googlePtms), qkToName('31'), logger); o(job.source.files).deepEquals([tif1, tif2]); o(cutTiffArgs.length).equals(0); @@ -91,7 +91,7 @@ o.spec('cog.vrt', () => { o('no cutline diff projection', async () => { job.projection = EpsgCode.Nztm2000; - const vrt = await CogVrt.buildVrt(tmpFolder, job, new Cutline(googleProj), qkToName('31'), logger); + const vrt = await CogVrt.buildVrt(tmpFolder, job, new Cutline(googlePtms), qkToName('31'), logger); o(job.source.files).deepEquals([tif1, tif2]); o(cutTiffArgs.length).equals(0); @@ -101,7 +101,7 @@ o.spec('cog.vrt', () => { }); o('fully within same projection', async () => { - const cutline = new Cutline(googleProj, await Cutline.loadCutline(testDir + '/kapiti.geojson'), -100); + const cutline = new Cutline(googlePtms, await Cutline.loadCutline(testDir + '/kapiti.geojson'), -100); const name = qkToName('311333222321113310'); @@ -115,7 +115,7 @@ o.spec('cog.vrt', () => { }); o('intersected cutline', async () => { - const cutline = new Cutline(googleProj, await Cutline.loadCutline(testDir + '/kapiti.geojson'), 20); + const cutline = new Cutline(googlePtms, await Cutline.loadCutline(testDir + '/kapiti.geojson'), 20); job.output.cutline = { blend: 20, source: 'cutline.json' }; const name = qkToName('311333222321113'); @@ -177,7 +177,7 @@ o.spec('cog.vrt', () => { return o; }); job.output.cutline = { blend: 10, source: 'cutline.json' }; - const cutline = new Cutline(googleProj, await Cutline.loadCutline(testDir + '/mana.geojson')); + const cutline = new Cutline(googlePtms, await Cutline.loadCutline(testDir + '/mana.geojson')); const vrt = await CogVrt.buildVrt(tmpFolder, job, cutline, qkToName('3131110001'), logger); diff --git a/packages/cli/src/cog/__test__/cutline.test.ts b/packages/cli/src/cog/__test__/cutline.test.ts index f76635546..1df79bb9b 100644 --- a/packages/cli/src/cog/__test__/cutline.test.ts +++ b/packages/cli/src/cog/__test__/cutline.test.ts @@ -7,16 +7,16 @@ import o from 'ospec'; import { Cutline } from '../cutline'; import { SourceMetadata } from '../types'; import { SourceTiffTestHelper } from './source.tiff.testhelper'; -import { qkToName } from '@basemaps/shared/build/tms/__test__/test.util'; +import { qkToName } from '@basemaps/shared/build/proj/__test__/test.util'; o.spec('cutline', () => { const testDir = `${__dirname}/../../../__test.assets__`; - const googleProj = ProjectionTileMatrixSet.get(EpsgCode.Google); - const nztmProj = ProjectionTileMatrixSet.get(EpsgCode.Nztm2000); + const googlePtms = ProjectionTileMatrixSet.get(EpsgCode.Google); + const nztmPtms = ProjectionTileMatrixSet.get(EpsgCode.Nztm2000); o.spec('filterSourcesForName', () => { o('fully within same projection', async () => { - const cutline = new Cutline(googleProj, await Cutline.loadCutline(testDir + '/kapiti.geojson'), -100); + const cutline = new Cutline(googlePtms, await Cutline.loadCutline(testDir + '/kapiti.geojson'), -100); const name = qkToName('311333222321113310'); @@ -34,7 +34,7 @@ o.spec('cutline', () => { }); o('loadCutline', async () => { - const cutline = new Cutline(googleProj, await Cutline.loadCutline(testDir + '/mana.geojson')); + const cutline = new Cutline(googlePtms, await Cutline.loadCutline(testDir + '/mana.geojson')); const geojson = cutline.toGeoJson(); const mp = geojson.features[0].geometry as MultiPolygon; const { coordinates } = mp; @@ -58,7 +58,7 @@ o.spec('cutline', () => { }); o('findCovering', async () => { - const cutline = new Cutline(googleProj, await Cutline.loadCutline(testDir + '/mana.geojson')); + const cutline = new Cutline(googlePtms, await Cutline.loadCutline(testDir + '/mana.geojson')); const bounds = SourceTiffTestHelper.tiffNztmBounds(); o(cutline.clipPoly.length).equals(2); @@ -72,7 +72,7 @@ o.spec('cutline', () => { const bounds = SourceTiffTestHelper.tiffNztmBounds(); o('nztm', () => { - const cutline = new Cutline(nztmProj); + const cutline = new Cutline(nztmPtms); const covering = cutline.optimizeCovering({ projection: EpsgCode.Nztm2000, @@ -122,7 +122,7 @@ o.spec('cutline', () => { }; const bounds = [tiff1, tiff2]; - const cutline = new Cutline(googleProj, await Cutline.loadCutline(testDir + '/kapiti.geojson'), 500); + const cutline = new Cutline(googlePtms, await Cutline.loadCutline(testDir + '/kapiti.geojson'), 500); const covering = cutline.optimizeCovering({ projection: EpsgCode.Nztm2000, @@ -158,7 +158,7 @@ o.spec('cutline', () => { }); o('boundary part inland, part coastal', async () => { - const cutline = new Cutline(googleProj, await Cutline.loadCutline(testDir + '/kapiti.geojson'), 500); + const cutline = new Cutline(googlePtms, await Cutline.loadCutline(testDir + '/kapiti.geojson'), 500); const bounds = [ { name: 'tiff1', @@ -197,7 +197,7 @@ o.spec('cutline', () => { }); o('low res', () => { - const cutline = new Cutline(googleProj); + const cutline = new Cutline(googlePtms); const covering = cutline.optimizeCovering({ projection: EpsgCode.Nztm2000, @@ -210,7 +210,7 @@ o.spec('cutline', () => { }); o('hi res', () => { - const covering2 = new Cutline(googleProj).optimizeCovering({ + const covering2 = new Cutline(googlePtms).optimizeCovering({ projection: EpsgCode.Nztm2000, bounds, resZoom: 18, diff --git a/packages/cli/src/cog/builder.ts b/packages/cli/src/cog/builder.ts index 449c17582..fbccf0e1a 100644 --- a/packages/cli/src/cog/builder.ts +++ b/packages/cli/src/cog/builder.ts @@ -43,7 +43,7 @@ export function guessProjection(wkt: string): Epsg | null { export class CogBuilder { q: pLimit.Limit; logger: LogType; - targetProj: ProjectionTileMatrixSet; + targetPtms: ProjectionTileMatrixSet; srcProj?: Epsg; // Prevent guessing spamming the logs @@ -52,8 +52,8 @@ export class CogBuilder { /** * @param concurrency number of requests to run at a time */ - constructor(targetProj: ProjectionTileMatrixSet, concurrency: number, logger: LogType, srcProj?: Epsg) { - this.targetProj = targetProj; + constructor(targetPtms: ProjectionTileMatrixSet, concurrency: number, logger: LogType, srcProj?: Epsg) { + this.targetPtms = targetPtms; this.logger = logger; this.q = pLimit(concurrency); this.srcProj = srcProj; @@ -136,7 +136,7 @@ export class CogBuilder { bands, bounds, pixelScale: resX, - resZoom: this.targetProj.getTiffResZoom(resX), + resZoom: this.targetPtms.getTiffResZoom(resX), }; } @@ -196,7 +196,7 @@ export class CogBuilder { CacheFolder, CacheVersion + createHash('sha256') - .update(this.targetProj.tms.projection.toString()) + .update(this.targetPtms.tms.projection.toString()) .update(tiffs.map((c) => c.name).join('\n')) .digest('hex'), ) + '.json'; @@ -204,7 +204,7 @@ export class CogBuilder { if (existsSync(cacheKey)) { this.logger.debug({ path: cacheKey }, 'MetadataCacheHit'); const metadata = (await FileOperatorSimple.readJson(cacheKey)) as SourceMetadata; - metadata.resZoom = this.targetProj.getTiffResZoom(metadata.pixelScale); + metadata.resZoom = this.targetPtms.getTiffResZoom(metadata.pixelScale); return metadata; } diff --git a/packages/cli/src/cog/cog.ts b/packages/cli/src/cog/cog.ts index 8bdec8391..fc2a336dd 100644 --- a/packages/cli/src/cog/cog.ts +++ b/packages/cli/src/cog/cog.ts @@ -49,20 +49,20 @@ export async function buildCogForName( const startTime = Date.now(); const { resZoom } = job.source; - const targetProj = ProjectionTileMatrixSet.get(job.projection); - const { tms } = targetProj; + const targetPtms = ProjectionTileMatrixSet.get(job.projection); + const { tms } = targetPtms; const tile = TileMatrixSet.nameToTile(name); const ul = tms.tileToSource(tile); const lr = tms.tileToSource({ x: tile.x + 1, y: tile.y + 1, z: tile.z }); - const blockSize = tms.tileSize * targetProj.blockFactor; - const alignmentLevels = targetProj.findAlignmentLevels(tile, job.source.pixelScale); + const blockSize = tms.tileSize * targetPtms.blockFactor; + const alignmentLevels = targetPtms.findAlignmentLevels(tile, job.source.pixelScale); const cogBuild = new GdalCogBuilder(vrtLocation, outputTiffPath, { bbox: [ul.x, ul.y, lr.x, lr.y], - projection: targetProj.tms.projection, + projection: targetPtms.tms.projection, tilingScheme: tilingScheme(job.projection), blockSize, targetRes: tms.pixelScale(resZoom), @@ -75,7 +75,7 @@ export async function buildCogForName( logger.info( { - imageSize: targetProj.getImagePixelWidth(tile, resZoom), + imageSize: targetPtms.getImagePixelWidth(tile, resZoom), name, tile, alignmentLevels, diff --git a/packages/cli/src/cog/cog.vrt.ts b/packages/cli/src/cog/cog.vrt.ts index 9cebfffeb..9430ad078 100644 --- a/packages/cli/src/cog/cog.vrt.ts +++ b/packages/cli/src/cog/cog.vrt.ts @@ -127,7 +127,7 @@ export const CogVrt = { gdalCommand.mount(tmpFolder); } - const tr = cutline.targetProj.tms.pixelScale(job.source.resZoom).toString(); + const tr = cutline.tms.pixelScale(job.source.resZoom).toString(); const defaultOps = ['-tr', tr, tr, '-tap']; diff --git a/packages/cli/src/cog/cutline.ts b/packages/cli/src/cog/cutline.ts index c041be1d6..579c4567d 100644 --- a/packages/cli/src/cog/cutline.ts +++ b/packages/cli/src/cog/cutline.ts @@ -5,6 +5,7 @@ import { CoveringFraction, MaxImagePixelWidth } from './constants'; import { CogJob, SourceMetadata } from './types'; import { clipMultipolygon, removeDegenerateEdges, polyContainsBounds } from './clipped.multipolygon'; import pc, { MultiPolygon, Ring, Polygon } from 'polygon-clipping'; +import { Projection } from '@basemaps/shared/build/proj/projection'; const { intersection, union } = pc; /** Padding to always apply to image boundies */ @@ -39,15 +40,15 @@ function addNonDupes(list: Tile[], addList: Tile[]): void { export class Cutline { clipPoly: MultiPolygon = []; - targetProj: ProjectionTileMatrixSet; + targetPtms: ProjectionTileMatrixSet; blend: number; - tms: TileMatrixSet; // convience to targetProj.tms + tms: TileMatrixSet; // convience to targetPtms.tms private srcPoly: MultiPolygon = []; /** * Create a Cutline instance from a `GeoJSON FeatureCollection`. - * @param targetProj the projection the COGs will be created in. + * @param targetPtms the projection the COGs will be created in. * @param clipPoly the optional cutline. The source imagery outline used by default. This * `FeatureCollection` is converted to one `MultiPolygon` with any holes removed and the @@ -55,9 +56,9 @@ export class Cutline { * @param blend How much blending to consider when working out boundaries. */ - constructor(targetProj: ProjectionTileMatrixSet, clipPoly?: FeatureCollection, blend = 0) { - this.targetProj = targetProj; - this.tms = targetProj.tms; + constructor(targetPtms: ProjectionTileMatrixSet, clipPoly?: FeatureCollection, blend = 0) { + this.targetPtms = targetPtms; + this.tms = targetPtms.tms; this.blend = blend; if (clipPoly == null) { return; @@ -65,7 +66,7 @@ export class Cutline { if (findGeoJsonProjection(clipPoly) !== Epsg.Wgs84) { throw new Error('Invalid geojson; CRS may not be set for cutline!'); } - const { fromWsg84 } = this.targetProj; + const { fromWsg84 } = this.targetPtms.proj; for (const { geometry } of clipPoly.features) { if (geometry.type === 'MultiPolygon') { for (const coords of geometry.coordinates) { @@ -96,16 +97,16 @@ export class Cutline { */ filterSourcesForName(name: string, job: CogJob): void { const tile = TileMatrixSet.nameToTile(name); - const { targetProj } = this; - const tileBounds = targetProj.tileToSourceBounds(tile); + const sourceCode = Projection.get(job.source.projection); + const targetCode = this.targetPtms.proj; + const tileBounds = this.targetPtms.tileToSourceBounds(tile); const tilePadded = this.padBounds(tileBounds, job.source.resZoom); let tileBoundsInSrcProj = tilePadded; - if (job.source.projection !== this.tms.projection.code) { + if (sourceCode !== targetCode) { // convert the padded quadKey to source projection ensuring fully enclosed - const sourceProj = ProjectionTileMatrixSet.get(job.source.projection); - const poly = targetProj.projectMultipolygon([tileBoundsInSrcProj.toPolygon()], sourceProj); + const poly = targetCode.projectMultipolygon([tileBoundsInSrcProj.toPolygon()], sourceCode); tileBoundsInSrcProj = Bounds.fromMultiPolygon(poly); } @@ -137,13 +138,13 @@ export class Cutline { const { resZoom } = sourceMetadata; // Look for the biggest tile size we are allowed to create. let minZ = resZoom - 1; - while (minZ > 1 && this.targetProj.getImagePixelWidth({ x: 0, y: 0, z: minZ }, resZoom) < MaxImagePixelWidth) { + while (minZ > 1 && this.targetPtms.getImagePixelWidth({ x: 0, y: 0, z: minZ }, resZoom) < MaxImagePixelWidth) { --minZ; } minZ = Math.max(1, minZ + 1); let tiles: Tile[] = []; - const { tms } = this.targetProj; + const { tms } = this.targetPtms; for (const tile of tms.topLevelTiles()) { // Don't make COGs with a tile.z < minZ. @@ -170,7 +171,7 @@ export class Cutline { * Convert JobCutline to geojson FeatureCollection */ toGeoJson(clipPoly = this.clipPoly): FeatureCollection { - const { toWsg84 } = this.targetProj; + const { toWsg84 } = this.targetPtms.proj; return GeoJson.toFeatureCollection([GeoJson.toFeatureMultiPolygon(clipPoly.map((p) => [p[0].map(toWsg84)]))]); } @@ -189,7 +190,7 @@ export class Cutline { minZ: number, coveringFraction: number, ): { tiles: Tile[]; fractionCovered: number } { - const clipBounds = this.targetProj.tileToSourceBounds(tile); + const clipBounds = this.targetPtms.tileToSourceBounds(tile); srcArea = clipMultipolygon(srcArea, clipBounds); @@ -242,8 +243,11 @@ export class Cutline { } // Convert polygon to target projection - const sourceProj = ProjectionTileMatrixSet.get(sourceMetadata.projection); - srcPoly = sourceProj.projectMultipolygon(srcPoly, this.targetProj) as MultiPolygon; + const sourceProj = Projection.get(sourceMetadata.projection); + const targetProj = this.targetPtms.proj; + if (sourceProj != targetProj) { + srcPoly = sourceProj.projectMultipolygon(srcPoly, targetProj) as MultiPolygon; + } this.srcPoly = srcPoly; if (this.clipPoly.length == 0) return; diff --git a/packages/cli/src/cog/job.ts b/packages/cli/src/cog/job.ts index b47035e87..7d450ca4c 100644 --- a/packages/cli/src/cog/job.ts +++ b/packages/cli/src/cog/job.ts @@ -20,6 +20,7 @@ import { getJobPath, makeTempFolder } from '../cli/folder'; import { GdalCogBuilderDefaults, GdalCogBuilderOptionsResampling } from '../gdal/gdal.config'; import { Cutline } from './cutline'; import { CogJob } from './types'; +import { Projection } from '@basemaps/shared/build/proj/projection'; export const MaxConcurrencyDefault = 50; @@ -123,9 +124,9 @@ export const CogJobFactory = { const builder = new CogBuilder(ctx.targetProjection, maxConcurrency, logger, ctx.override?.projection); const metadata = await builder.build(tiffSource, cutline); - const sourceProjection = ProjectionTileMatrixSet.get(metadata.projection); + const sourceProjection = Projection.get(metadata.projection); - const { tms } = cutline.targetProj; + const { tms } = cutline.targetPtms; const files = metadata.files.sort(Bounds.compareArea); if (files.length > 0) { @@ -212,7 +213,11 @@ export const CogJobFactory = { await outputFs.writeJson(geoJsonSourceOutput, sourceProjection.toGeoJson(metadata.bounds), logger); const geoJsonCoveringOutput = getJobPath(job, `covering.geojson`); - await outputFs.writeJson(geoJsonCoveringOutput, ctx.targetProjection.toGeoJson(metadata.files), logger); + await outputFs.writeJson( + geoJsonCoveringOutput, + ctx.targetProjection.proj.toGeoJson(metadata.files), + logger, + ); if (ctx.batch) { await ActionBatchJob.batchJob(jobFile, true, undefined, logger); diff --git a/packages/shared/src/aws/__test__/tile.metadata.table.test.ts b/packages/shared/src/aws/__test__/tile.metadata.table.test.ts index 0c878605b..cdd117570 100644 --- a/packages/shared/src/aws/__test__/tile.metadata.table.test.ts +++ b/packages/shared/src/aws/__test__/tile.metadata.table.test.ts @@ -3,7 +3,7 @@ import { round } from '@basemaps/test/build/rounding'; import * as AWS from 'aws-sdk'; import o from 'ospec'; import { Const } from '../../const'; -import { qkToNamedBounds } from '../../tms/__test__/test.util'; +import { qkToNamedBounds } from '../../proj/__test__/test.util'; import { TileMetadataTable } from '../tile.metadata'; import { TileMetadataImageRule, TileMetadataImageryRecordV1, TileMetadataSetRecord } from '../tile.metadata.base'; diff --git a/packages/shared/src/index.ts b/packages/shared/src/index.ts index 17ee5c743..0b0f34fd5 100644 --- a/packages/shared/src/index.ts +++ b/packages/shared/src/index.ts @@ -9,7 +9,7 @@ export { V, VNode, VNodeElement, VNodeText } from './vdom'; export { VNodeParser } from './vdom.parse'; export { CompositeError } from './composite.error'; export { LoggerFatalError } from './logger.fatal.error'; -export * from './tms/projection.tile.matrix.set'; +export * from './proj/projection.tile.matrix.set'; export * from './aws/tile.metadata.base'; export * from './aws/tile.metadata'; diff --git a/packages/shared/src/proj/__test__/projection.test.ts b/packages/shared/src/proj/__test__/projection.test.ts new file mode 100644 index 000000000..f678cef71 --- /dev/null +++ b/packages/shared/src/proj/__test__/projection.test.ts @@ -0,0 +1,68 @@ +import { Bounds, EpsgCode } from '@basemaps/geo'; +import { round } from '@basemaps/test/build/rounding'; +import o from 'ospec'; +import { Projection } from '../projection'; +import { qkToNamedBounds } from './test.util'; + +o.spec('Projection', () => { + const googleProj = Projection.get(EpsgCode.Google); + const nztmProj = Projection.get(EpsgCode.Nztm2000); + + o('should convert to 2193', () => { + if (nztmProj == null) { + throw new Error('Failed to init proj:2193'); + } + const output = nztmProj.toWsg84([1180000, 4758000]); + o(round(output, 6)).deepEquals([167.454458, -47.197075]); + + const reverse = nztmProj.fromWsg84(output); + o(round(reverse, 2)).deepEquals([1180000, 4758000]); + }); + + o('toGeoJson', () => { + const geojson = googleProj.toGeoJson(qkToNamedBounds(['31', '33'])); + + const { features } = geojson; + + o(features.length).equals(2); + + o(features[0].properties).deepEquals({ name: '2-3-2' }); + o(features[1].properties).deepEquals({ name: '2-3-3' }); + const { geometry } = features[0]!; + const coords = geometry.type === 'Polygon' ? geometry.coordinates : null; + o(round(coords![0], 8)).deepEquals([ + [90, -66.51326044], + [90, 0], + [180, 0], + [180, -66.51326044], + [90, -66.51326044], + ]); + }); + + o('projectMultipolygon', () => { + const poly = [ + Bounds.fromBbox([ + 18494091.86765497, + -6051366.655280836, + 19986142.659781612, + -4016307.214216303, + ]).toPolygon(), + ]; + + o(googleProj.projectMultipolygon(poly, googleProj)).equals(poly); + + const ans = googleProj.projectMultipolygon(poly, nztmProj); + + o(round(ans, 4)).deepEquals([ + [ + [ + [1084733.8967, 4698018.9435], + [964788.1197, 6226878.2808], + [2204979.5633, 6228860.047], + [2090794.171, 4700144.6365], + [1084733.8967, 4698018.9435], + ], + ], + ]); + }); +}); diff --git a/packages/shared/src/tms/__test__/projection.tile.matrix.set.test.ts b/packages/shared/src/proj/__test__/projection.tile.matrix.set.test.ts similarity index 71% rename from packages/shared/src/tms/__test__/projection.tile.matrix.set.test.ts rename to packages/shared/src/proj/__test__/projection.tile.matrix.set.test.ts index 93cd60e64..b3c7a1b2c 100644 --- a/packages/shared/src/tms/__test__/projection.tile.matrix.set.test.ts +++ b/packages/shared/src/proj/__test__/projection.tile.matrix.set.test.ts @@ -4,7 +4,6 @@ import { Approx } from '@basemaps/test'; import { round } from '@basemaps/test/build/rounding'; import o from 'ospec'; import { ProjectionTileMatrixSet } from '../projection.tile.matrix.set'; -import { qkToNamedBounds } from './test.util'; const TileSize = 256; @@ -26,53 +25,42 @@ function getPixelsFromTile(x: number, y: number): Bounds { } o.spec('ProjectionTileMatrixSet', () => { - const googleProj = ProjectionTileMatrixSet.get(EpsgCode.Google); - const nztmProj = ProjectionTileMatrixSet.get(EpsgCode.Nztm2000); + const googlePtms = ProjectionTileMatrixSet.get(EpsgCode.Google); + const nztmPtms = ProjectionTileMatrixSet.get(EpsgCode.Nztm2000); o('getTiffResZoom', () => { - o(googleProj.getTiffResZoom(10)).equals(14); - o(googleProj.getTiffResZoom(10, 2)).equals(13); - o(googleProj.getTiffResZoom(0.075)).equals(21); + o(googlePtms.getTiffResZoom(10)).equals(14); + o(googlePtms.getTiffResZoom(10, 2)).equals(13); + o(googlePtms.getTiffResZoom(0.075)).equals(21); - o(nztmProj.getTiffResZoom(10)).equals(10); - o(nztmProj.getTiffResZoom(10, 2)).equals(9); - o(nztmProj.getTiffResZoom(0.075)).equals(16); + o(nztmPtms.getTiffResZoom(10)).equals(10); + o(nztmPtms.getTiffResZoom(10, 2)).equals(9); + o(nztmPtms.getTiffResZoom(0.075)).equals(16); }); o('getTileSize', async () => { - o(googleProj.getImagePixelWidth({ x: 0, y: 0, z: 5 }, 10)).equals(16384); - o(googleProj.getImagePixelWidth({ x: 0, y: 0, z: 13 }, 20)).equals(65536); + o(googlePtms.getImagePixelWidth({ x: 0, y: 0, z: 5 }, 10)).equals(16384); + o(googlePtms.getImagePixelWidth({ x: 0, y: 0, z: 13 }, 20)).equals(65536); - o(nztmProj.getImagePixelWidth({ x: 0, y: 0, z: 5 }, 10)).equals(20480); - o(nztmProj.getImagePixelWidth({ x: 0, y: 0, z: 13 }, 16)).equals(5120); + o(nztmPtms.getImagePixelWidth({ x: 0, y: 0, z: 5 }, 10)).equals(20480); + o(nztmPtms.getImagePixelWidth({ x: 0, y: 0, z: 13 }, 16)).equals(5120); }); o('findAlignmentLevels', () => { - o(googleProj.findAlignmentLevels({ x: 2, y: 0, z: 5 }, 0.075)).equals(15); - o(googleProj.findAlignmentLevels({ x: 2, y: 0, z: 5 }, 0.5)).equals(13); - o(googleProj.findAlignmentLevels({ x: 2, y: 0, z: 3 }, 1)).equals(14); - o(googleProj.findAlignmentLevels({ x: 2, y: 0, z: 8 }, 10)).equals(5); - o(googleProj.findAlignmentLevels({ x: 2, y: 0, z: 14 }, 10)).equals(0); - - o(nztmProj.findAlignmentLevels({ x: 2, y: 0, z: 1 }, 0.075)).equals(14); - o(nztmProj.findAlignmentLevels({ x: 2, y: 0, z: 5 }, 0.5)).equals(8); - o(nztmProj.findAlignmentLevels({ x: 2, y: 0, z: 3 }, 7)).equals(6); - o(nztmProj.findAlignmentLevels({ x: 2, y: 0, z: 8 }, 14)).equals(0); - }); - - o('should convert to 2193', () => { - if (nztmProj == null) { - throw new Error('Failed to init proj:2193'); - } - const output = nztmProj.toWsg84([1180000, 4758000]); - o(round(output, 6)).deepEquals([167.454458, -47.197075]); - - const reverse = nztmProj.fromWsg84(output); - o(round(reverse, 2)).deepEquals([1180000, 4758000]); + o(googlePtms.findAlignmentLevels({ x: 2, y: 0, z: 5 }, 0.075)).equals(15); + o(googlePtms.findAlignmentLevels({ x: 2, y: 0, z: 5 }, 0.5)).equals(13); + o(googlePtms.findAlignmentLevels({ x: 2, y: 0, z: 3 }, 1)).equals(14); + o(googlePtms.findAlignmentLevels({ x: 2, y: 0, z: 8 }, 10)).equals(5); + o(googlePtms.findAlignmentLevels({ x: 2, y: 0, z: 14 }, 10)).equals(0); + + o(nztmPtms.findAlignmentLevels({ x: 2, y: 0, z: 1 }, 0.075)).equals(14); + o(nztmPtms.findAlignmentLevels({ x: 2, y: 0, z: 5 }, 0.5)).equals(8); + o(nztmPtms.findAlignmentLevels({ x: 2, y: 0, z: 3 }, 7)).equals(6); + o(nztmPtms.findAlignmentLevels({ x: 2, y: 0, z: 8 }, 14)).equals(0); }); o('tileToSourceBounds', () => { - o(round(googleProj.tileToSourceBounds(QuadKey.toTile('3120123')).toJson(), 8)).deepEquals({ + o(round(googlePtms.tileToSourceBounds(QuadKey.toTile('3120123')).toJson(), 8)).deepEquals({ x: 11584184.51067502, y: -6261721.35712163, width: 313086.06785608, @@ -81,59 +69,12 @@ o.spec('ProjectionTileMatrixSet', () => { }); o('tileCenterToLatLon', () => { - o(round(googleProj.tileCenterToLatLon(QuadKey.toTile('3120123')), 8)).deepEquals({ + o(round(googlePtms.tileCenterToLatLon(QuadKey.toTile('3120123')), 8)).deepEquals({ lat: -47.98992167, lon: 105.46875, }); }); - o('toGeoJson', () => { - const geojson = googleProj.toGeoJson(qkToNamedBounds(['31', '33'], googleProj)); - - const { features } = geojson; - - o(features.length).equals(2); - - o(features[0].properties).deepEquals({ name: '2-3-2' }); - o(features[1].properties).deepEquals({ name: '2-3-3' }); - const { geometry } = features[0]!; - const coords = geometry.type === 'Polygon' ? geometry.coordinates : null; - o(round(coords![0], 8)).deepEquals([ - [90, -66.51326044], - [90, 0], - [180, 0], - [180, -66.51326044], - [90, -66.51326044], - ]); - }); - - o('projectMultipolygon', () => { - const poly = [ - Bounds.fromBbox([ - 18494091.86765497, - -6051366.655280836, - 19986142.659781612, - -4016307.214216303, - ]).toPolygon(), - ]; - - o(googleProj.projectMultipolygon(poly, googleProj)).equals(poly); - - const ans = googleProj.projectMultipolygon(poly, nztmProj); - - o(round(ans, 4)).deepEquals([ - [ - [ - [1084733.8967, 4698018.9435], - [964788.1197, 6226878.2808], - [2204979.5633, 6228860.047], - [2090794.171, 4700144.6365], - [1084733.8967, 4698018.9435], - ], - ], - ]); - }); - o.spec('TilingBounds', () => { // Approximate bounding box of new zealand const tifBoundingBox: [number, number, number, number] = [ diff --git a/packages/shared/src/tms/__test__/test.util.ts b/packages/shared/src/proj/__test__/test.util.ts similarity index 85% rename from packages/shared/src/tms/__test__/test.util.ts rename to packages/shared/src/proj/__test__/test.util.ts index 683b4d06a..4e123df2c 100644 --- a/packages/shared/src/tms/__test__/test.util.ts +++ b/packages/shared/src/proj/__test__/test.util.ts @@ -8,8 +8,8 @@ export function qkToName(qk: string): string { export function qkToNamedBounds( quadKeys: string[], - proj = ProjectionTileMatrixSet.get(EpsgCode.Google), + ptms = ProjectionTileMatrixSet.get(EpsgCode.Google), ): NamedBounds[] { - const { tms } = proj; + const { tms } = ptms; return quadKeys.map((qk) => Object.assign({ name: qkToName(qk) }, tms.tileToSourceBounds(QuadKey.toTile(qk)))); } diff --git a/packages/shared/src/tms/citm2000.ts b/packages/shared/src/proj/citm2000.ts similarity index 100% rename from packages/shared/src/tms/citm2000.ts rename to packages/shared/src/proj/citm2000.ts diff --git a/packages/shared/src/tms/nztm2000.ts b/packages/shared/src/proj/nztm2000.ts similarity index 100% rename from packages/shared/src/tms/nztm2000.ts rename to packages/shared/src/proj/nztm2000.ts diff --git a/packages/shared/src/tms/projection.tile.matrix.set.ts b/packages/shared/src/proj/projection.tile.matrix.set.ts similarity index 57% rename from packages/shared/src/tms/projection.tile.matrix.set.ts rename to packages/shared/src/proj/projection.tile.matrix.set.ts index d9f34d022..31d5bb1d4 100644 --- a/packages/shared/src/tms/projection.tile.matrix.set.ts +++ b/packages/shared/src/proj/projection.tile.matrix.set.ts @@ -1,46 +1,29 @@ -import { Bounds, Epsg, EpsgCode, GeoJson, Tile, TileMatrixSet } from '@basemaps/geo'; +import { Bounds, EpsgCode, Tile, TileMatrixSet } from '@basemaps/geo'; import { GoogleTms } from '@basemaps/geo/build/tms/google'; import { Nztm2000Tms } from '@basemaps/geo/build/tms/nztm2000'; -import { Position } from 'geojson'; -import Proj from 'proj4'; -import { NamedBounds } from '../aws/tile.metadata.base'; -import { CompositeError } from '../composite.error'; -import { Citm2000 } from './citm2000'; -import { Nztm2000 } from './nztm2000'; +import { Projection } from './projection'; export interface LatLon { lat: number; lon: number; } -Proj.defs(Epsg.Nztm2000.toEpsgString(), Nztm2000); -Proj.defs(Epsg.Citm2000.toEpsgString(), Citm2000); - const CodeMap = new Map(); export class ProjectionTileMatrixSet { /** The underlying TileMatrixSet */ - tms: TileMatrixSet; + public readonly tms: TileMatrixSet; + public readonly proj: Projection; /** Used to calculate `BlockSize = blockFactor * tms.tileSize` for generating COGs */ blockFactor: number; - /** Transform coordinates to and from Wsg84 */ - private projection: proj4.Converter; - /** * Wrapper around TileMatrixSet with utilities for converting Points and Polygons */ constructor(tms: TileMatrixSet, blockFactor = 2) { this.tms = tms; this.blockFactor = blockFactor; - try { - this.projection = Proj(tms.projection.toEpsgString(), Epsg.Wgs84.toEpsgString()); - } catch (err) { - throw new CompositeError( - `Failed to create projection: ${tms.projection.toEpsgString()}, ${Epsg.Wgs84.toEpsgString()}`, - err, - ); - } + this.proj = Projection.get(tms.projection.code); } /** @@ -64,35 +47,6 @@ export class ProjectionTileMatrixSet { return (epsgCode && CodeMap.get(epsgCode)) ?? null; } - /** - * Project the points in a MultiPolygon array to the `targetProjection`. - - * @return if multipoly is not projected return it verbatim otherwise creates a new multi - * polygon - */ - projectMultipolygon(multipoly: Position[][][], targetProjection: ProjectionTileMatrixSet): Position[][][] { - if (targetProjection.tms.projection.code === this.tms.projection.code) return multipoly; - - const { toWsg84 } = this; - const { fromWsg84 } = targetProjection; - - return multipoly.map((poly) => poly.map((ring) => ring.map((p) => fromWsg84(toWsg84(p))))); - } - - /** - * Convert source `[x, y]` coordinates to `[lon, lat]` - */ - get toWsg84(): (coordinates: Position) => Position { - return this.projection.forward; - } - - /** - * Convert `[lon, lat]` coordinates to source `[x, y]` - */ - get fromWsg84(): (coordinates: Position) => Position { - return this.projection.inverse; - } - /** * Find the closest zoom level to `resX` (pixels per meter) that is at least as good as `resX`. @@ -124,31 +78,10 @@ export class ProjectionTileMatrixSet { */ tileCenterToLatLon(tile: Tile): LatLon { const point = this.tms.tileToSource({ x: tile.x + 0.5, y: tile.y + 0.5, z: tile.z }); - const [lon, lat] = this.toWsg84([point.x, point.y]); + const [lon, lat] = this.proj.toWsg84([point.x, point.y]); return { lat, lon }; } - /** Convert a tile covering to a GeoJSON FeatureCollection */ - toGeoJson(files: NamedBounds[]): GeoJSON.FeatureCollection { - const { toWsg84 } = this; - const polygons: GeoJSON.Feature[] = []; - for (const bounds of files) { - const polygon = GeoJson.toFeaturePolygon( - [ - Bounds.fromJson(bounds) - .toPolygon()[0] - .map((p) => toWsg84(p)), - ], - { - name: bounds.name, - }, - ); - polygons.push(polygon); - } - - return GeoJson.toFeatureCollection(polygons); - } - /** * Find the number of alignment levels required to render the tile. Min 1 * diff --git a/packages/shared/src/proj/projection.ts b/packages/shared/src/proj/projection.ts new file mode 100644 index 000000000..8d8f9ae42 --- /dev/null +++ b/packages/shared/src/proj/projection.ts @@ -0,0 +1,116 @@ +import { Bounds, Epsg, EpsgCode, GeoJson } from '@basemaps/geo'; +import { Position } from 'geojson'; +import Proj from 'proj4'; +import { NamedBounds } from '../aws/tile.metadata.base'; +import { CompositeError } from '../composite.error'; +import { Citm2000 } from './citm2000'; +import { Nztm2000 } from './nztm2000'; + +Proj.defs(Epsg.Nztm2000.toEpsgString(), Nztm2000); +Proj.defs(Epsg.Citm2000.toEpsgString(), Citm2000); + +const CodeMap = new Map(); + +export class Projection { + epsg: Epsg; + + /** Transform coordinates to and from Wsg84 */ + private projection: proj4.Converter; + + /** + * Wrapper around TileMatrixSet with utilities for converting Points and Polygons + */ + private constructor(epsg: Epsg) { + this.epsg = epsg; + try { + this.projection = Proj(epsg.toEpsgString(), Epsg.Wgs84.toEpsgString()); + } catch (err) { + throw new CompositeError( + `Failed to create projection: ${epsg.toEpsgString()}, ${Epsg.Wgs84.toEpsgString()}`, + err, + ); + } + } + + /** + * Get the Projection instance for a specified code, + * + * throws a exception if the code is not recognized + * + * @param epsgCode + */ + static get(epsgCode: EpsgCode): Projection { + let proj = CodeMap.get(epsgCode); + if (proj != null) return proj; + const epsg = Epsg.tryGet(epsgCode); + if (epsg == null) { + throw new Error(`Invalid projection: ${epsgCode}`); + } + proj = new Projection(epsg); + CodeMap.set(epsgCode, proj); + return proj; + } + + /** + * Try to find a corresponding Projection for a number + * @param epsgCode + */ + static tryGet(epsgCode?: EpsgCode): Projection | null { + if (epsgCode == null) return null; + try { + return this.get(epsgCode); + } catch (err) { + return null; + } + } + + /** + * Project the points in a MultiPolygon array to the `targetProjection`. + + * @return if multipoly is not projected return it verbatim otherwise creates a new multi + * polygon + */ + projectMultipolygon(multipoly: Position[][][], targetProjection: Projection): Position[][][] { + if (targetProjection.epsg.code === this.epsg.code) return multipoly; + + const { toWsg84 } = this; + const { fromWsg84 } = targetProjection; + + return multipoly.map((poly) => poly.map((ring) => ring.map((p) => fromWsg84(toWsg84(p))))); + } + + /** + * Convert source `[x, y]` coordinates to `[lon, lat]` + */ + get toWsg84(): (coordinates: Position) => Position { + return this.projection.forward; + } + + /** + * Convert `[lon, lat]` coordinates to source `[x, y]` + */ + get fromWsg84(): (coordinates: Position) => Position { + return this.projection.inverse; + } + + /** Convert a tile covering to a GeoJSON FeatureCollection */ + toGeoJson(files: NamedBounds[]): GeoJSON.FeatureCollection { + const { toWsg84 } = this; + const polygons: GeoJSON.Feature[] = []; + for (const bounds of files) { + const polygon = GeoJson.toFeaturePolygon( + [ + Bounds.fromJson(bounds) + .toPolygon()[0] + .map((p) => toWsg84(p)), + ], + { + name: bounds.name, + }, + ); + polygons.push(polygon); + } + + return GeoJson.toFeatureCollection(polygons); + } +}