From cd33f1611826f42b7eeae65aa084d19d1fdaf767 Mon Sep 17 00:00:00 2001 From: Geoff Jacobsen Date: Sat, 27 Jun 2020 17:39:09 +1200 Subject: [PATCH] fix(cli): mitigate polygon intersection errors Use clipline to reduce polygon when we don't care about degenerate edges. When we do care call intersection to remove the degenerate edges. Store image bounds not polygons in SourceMetadata and also store in CogJob. This reduces floating-point errors when combining source imagery bounds into composite polygons. Note: The Martinez polygon-clipping libraries can suffer from errors when lines align. Other intersection libraries are intractable for the size of coastline polygons. A solution is to call intersection after first calling clipline. --- packages/cli/package.json | 1 + packages/cli/src/cli/cogify/action.cog.ts | 4 +- packages/cli/src/cog/__test__/builder.test.ts | 99 ++++-- .../cog/__test__/clipped.multipolygon.test.ts | 299 ++++++++++++++++++ packages/cli/src/cog/__test__/cog.vrt.test.ts | 125 +++++--- packages/cli/src/cog/__test__/cutline.test.ts | 131 ++++++-- .../cog/__test__/source.tiff.testhelper.ts | 48 ++- packages/cli/src/cog/builder.ts | 44 +-- packages/cli/src/cog/clipped.multipolygon.ts | 31 ++ packages/cli/src/cog/cog.vrt.ts | 27 +- packages/cli/src/cog/cutline.ts | 151 ++++----- packages/cli/src/cog/job.ts | 8 +- packages/cli/src/cog/types.ts | 4 +- packages/cli/src/lineclip.d.ts | 3 + packages/geo/src/__tests__/bounds.test.ts | 30 +- packages/geo/src/bounds.ts | 39 +++ packages/geo/src/geo.json.ts | 12 - .../projection.tile.matrix.set.test.ts | 25 ++ .../src/tms/projection.tile.matrix.set.ts | 38 +-- yarn.lock | 5 + 20 files changed, 829 insertions(+), 295 deletions(-) create mode 100644 packages/cli/src/cog/__test__/clipped.multipolygon.test.ts create mode 100644 packages/cli/src/cog/clipped.multipolygon.ts create mode 100644 packages/cli/src/lineclip.d.ts diff --git a/packages/cli/package.json b/packages/cli/package.json index 6fa9474e3e..a3b6d0ad9d 100644 --- a/packages/cli/package.json +++ b/packages/cli/package.json @@ -25,6 +25,7 @@ "@cogeotiff/source-url": "^2.1.1", "@rushstack/ts-command-line": "^4.3.13", "ansi-colors": "^4.1.1", + "lineclip": "^1.1.5", "p-limit": "^3.0.1", "polygon-clipping": "^0.14.3", "pretty-json-log": "^0.3.1" diff --git a/packages/cli/src/cli/cogify/action.cog.ts b/packages/cli/src/cli/cogify/action.cog.ts index 471925df89..0a0b1ad10f 100644 --- a/packages/cli/src/cli/cogify/action.cog.ts +++ b/packages/cli/src/cli/cogify/action.cog.ts @@ -101,8 +101,6 @@ export class ActionCogCreate extends CommandLineAction { const tmpFolder = await makeTempFolder(`basemaps-${job.id}-${CliId}`); try { - const sourceGeo = (await inFp.readJson(getJobPath(job, 'source.geojson'))) as FeatureCollection; - let cutlineJson: FeatureCollection | undefined; if (job.output.cutline != null) { const cutlinePath = getJobPath(job, 'cutline.geojson.gz'); @@ -117,7 +115,7 @@ export class ActionCogCreate extends CommandLineAction { job.output.cutline?.blend, ); - const tmpVrtPath = await CogVrt.buildVrt(tmpFolder, job, sourceGeo, cutline, name, logger); + const tmpVrtPath = await CogVrt.buildVrt(tmpFolder, job, cutline, name, logger); if (tmpVrtPath == null) { logger.warn({ name }, 'NoMatchingSourceImagery'); diff --git a/packages/cli/src/cog/__test__/builder.test.ts b/packages/cli/src/cog/__test__/builder.test.ts index c21f624838..96d433df8b 100644 --- a/packages/cli/src/cog/__test__/builder.test.ts +++ b/packages/cli/src/cog/__test__/builder.test.ts @@ -1,7 +1,8 @@ -import { Epsg, EpsgCode } from '@basemaps/geo'; +import { Bounds, Epsg, EpsgCode } from '@basemaps/geo'; import { LogConfig, ProjectionTileMatrixSet } from '@basemaps/shared'; -import { round } from '@basemaps/test/build/rounding'; -import { CogTiff, TiffTagGeo } from '@cogeotiff/core'; +import { CogTiff } from '@cogeotiff/core'; +import { CogSourceAwsS3 } from '@cogeotiff/source-aws'; +import { CogSourceFile } from '@cogeotiff/source-file'; import o from 'ospec'; import { CogBuilder, guessProjection } from '../builder'; @@ -28,35 +29,71 @@ o.spec('Builder', () => { o.spec('tiff', () => { const googleBuilder = new CogBuilder(ProjectionTileMatrixSet.get(EpsgCode.Google), 1, LogConfig.get()); - const tiff = { - source: { name: 'test1.tiff' }, - getImage(n: number): any { - if (n != 0) return null; - return { - bbox: [1492000, 6198000, 1492000 + 24000, 6198000 + 36000], - valueGeo(key: number): any { - if (key === TiffTagGeo.ProjectedCSTypeGeoKey) return EpsgCode.Nztm2000; + const origInit = CogTiff.prototype.init; + const origGetImage = CogTiff.prototype.getImage; + + o.after(() => { + CogTiff.prototype.init = origInit; + CogTiff.prototype.getImage = origGetImage; + }); + + o('bounds', async () => { + const localTiff = new CogSourceFile('local/file.tiff'); + const s3Tiff = new CogSourceAwsS3('bucket', 's3://file.tiff'); + + const imageLocal = { + resolution: [0.1], + value: (): any => [1], + valueGeo: (): any => EpsgCode.Nztm2000, + bbox: Bounds.fromJson({ + x: 1492000, + y: 6198000, + width: 24000, + height: 36000, + }).toBbox(), + }; + + const imageS3 = { + resolution: [0.1], + value: (): any => [1], + valueGeo: (): any => EpsgCode.Nztm2000, + bbox: Bounds.fromJson({ + x: 1492000 + 24000, + y: 6198000, + width: 24000, + height: 36000, + }).toBbox(), + }; + + CogTiff.prototype.init = o.spy() as any; + CogTiff.prototype.getImage = function (): any { + return this.source == localTiff ? imageLocal : imageS3; + }; + + const ans = await googleBuilder.bounds([localTiff, s3Tiff]); + + o(ans).deepEquals({ + projection: 2193, + nodata: 1, + bands: 1, + bounds: [ + { + x: 1492000, + y: 6198000, + width: 24000, + height: 36000, + name: 'local/file.tiff', + }, + { + x: 1516000, + y: 6198000, + width: 24000, + height: 36000, + name: 's3://bucket/s3://file.tiff', }, - }; - }, - } as CogTiff; - - o('getTifBounds', () => { - o(round(googleBuilder.getTifBounds(tiff), 2)).deepEquals({ - type: 'Feature', - geometry: { - type: 'Polygon', - coordinates: [ - [ - [171.83, -34.03], - [171.83, -34.35], - [172.09, -34.36], - [172.09, -34.03], - [171.83, -34.03], - ], - ], - }, - properties: { tiff: 'test1.tiff' }, + ], + pixelScale: 0.1, + resZoom: 21, }); }); }); diff --git a/packages/cli/src/cog/__test__/clipped.multipolygon.test.ts b/packages/cli/src/cog/__test__/clipped.multipolygon.test.ts new file mode 100644 index 0000000000..5e947e6ba4 --- /dev/null +++ b/packages/cli/src/cog/__test__/clipped.multipolygon.test.ts @@ -0,0 +1,299 @@ +import { Bounds, GeoJson } from '@basemaps/geo'; +import o from 'ospec'; +import { MultiPolygon } from 'polygon-clipping'; +import { clipMultipolygon, removeDegenerateEdges, polyContainsBounds } from '../clipped.multipolygon'; +import { round } from '@basemaps/test/build/rounding'; + +export function writeGeoJson(fn: string, poly: MultiPolygon): void { + require('fs').writeFileSync(fn, JSON.stringify(GeoJson.toFeatureCollection([GeoJson.toFeatureMultiPolygon(poly)]))); +} + +o.spec('clipped.multipolygon', () => { + const polys: MultiPolygon = [ + [ + [ + [-4, 2], + [-2, 1], + [-4, -1], + [1, -5], + [2, -5], + [2, -2], + [4, -5], + [6, -2], + [2, 0], + [6, 3], + [2, 7], + [0, 3], + [-4, 6], + [-4, 10], + [10, 10], + [10, -10], + [-10, -10], + [-4, 2], + ], + ], + [ + [ + [2, 3], + [3, 3], + [3, 5], + [2, 5], + [2, 3], + ], + ], + ]; + + o('polyContainsBounds', () => { + o(polyContainsBounds(polys, Bounds.fromBbox([-6, -6, -3, -3]))).equals(true); + + o(polyContainsBounds(polys, Bounds.fromBbox([-3, -4, 4, 4]))).equals(false); + o(polyContainsBounds(polys, new Bounds(6, -8, 5, 5))).equals(false); + }); + + o('clipMultipolygon', () => { + const bbox = Bounds.fromBbox([-3, -4, 4, 4]); + + const cp = clipMultipolygon(polys, bbox); + + o(cp).deepEquals([ + [ + [ + [-3, -4], + [-3, 1.5], + [-2, 1], + [-3, 0], + [-3, -1.8], + [-0.25, -4], + [2, -4], + [2, -2], + [3.333333333333333, -4], + [4, -4], + [4, -1], + [2, 0], + [4, 1.5], + [4, 4], + [0.5, 4], + [0, 3], + [-1.3333333333333333, 4], + [4, 4], + [4, -4], + [-3, -4], + ], + ], + [ + [ + [2, 3], + [3, 3], + [3, 4], + [2, 4], + [2, 3], + ], + ], + ]); + }); + + o.spec('removeDegenerateEdges', () => { + o('simple', () => { + const bbox = Bounds.fromBbox([-3, 1, 3, 3]); + + const cp: MultiPolygon = [ + [ + [ + [-2, 2], + [-1.5, 1], + [-0.5, 1], + [0, 2], + [0.5, 1], + [-2, 1], + [-2, 2], + ], + ], + ]; + + const ans = removeDegenerateEdges(cp, bbox); + + o(ans).deepEquals([ + [ + [ + [-2, 1], + [-1.5, 1], + [-2, 2], + [-2, 1], + ], + ], + [ + [ + [-0.5, 1], + [0.5, 1], + [0, 2], + [-0.5, 1], + ], + ], + ]); + }); + + o('three bumps', () => { + const degen: MultiPolygon = [ + [ + [ + [0, 0], + [159, 0], + [462, -3547], + [840, -3547], + [838, -3500], + [984, -3500], + [992, -3547], + [999, -3547], + [1013, -3444], + [1026, -3441], + [1045, -3442], + [1050, -3547], + [-693, -3547], + [-693, -3000], + [-2367, -3000], + [-2372, -3547], + [-2513, -3547], + [-2528, 0], + [0, 0], + ], + ], + ]; + + const bbox = Bounds.fromBbox([-3000, -3547, 3000, 1000]); + + const ans = removeDegenerateEdges(degen, bbox); + + o(ans).deepEquals([ + [ + [ + [-2528, 0], + [-2513, -3547], + [-2372, -3547], + [-2367, -3000], + [-693, -3000], + [-693, -3547], + [462, -3547], + [159, 0], + [-2528, 0], + ], + ], + [ + [ + [838, -3500], + [840, -3547], + [992, -3547], + [984, -3500], + [838, -3500], + ], + ], + [ + [ + [999, -3547], + [1050, -3547], + [1045, -3442], + [1026, -3441], + [1013, -3444], + [999, -3547], + ], + ], + ]); + }); + + o('loop', () => { + const bbox = Bounds.fromBbox([0, 3, 10, 8]); + + const orig: MultiPolygon = [ + [ + [ + [8, 6], + [2, 6], + [2, 1], + [4, 1], + [4, 4], + [6, 4], + [6, 1], + [8, 1], + [8, 6], + ], + ], + ]; + + const cp = clipMultipolygon(orig, bbox); + const ans = removeDegenerateEdges(cp, bbox); + + o(ans).deepEquals([ + [ + [ + [2, 3], + [4, 3], + [4, 4], + [6, 4], + [6, 3], + [8, 3], + [8, 6], + [2, 6], + [2, 3], + ], + ], + ]); + }); + + o('complex', () => { + const bbox = Bounds.fromBbox([-3, -4, 4, 4]); + const cp = clipMultipolygon(polys, bbox); + const ans = removeDegenerateEdges(cp, bbox); + + o(round(ans, 2)).deepEquals([ + [ + [ + [-3, -4], + [-0.25, -4], + [-3, -1.8], + [-3, -4], + ], + ], + [ + [ + [-3, 0], + [-2, 1], + [-3, 1.5], + [-3, 0], + ], + ], + [ + [ + [-1.33, 4], + [0, 3], + [0.5, 4], + [-1.33, 4], + ], + ], + [ + [ + [2, -4], + [3.33, -4], + [2, -2], + [2, -4], + ], + ], + [ + [ + [2, 0], + [4, -1], + [4, 1.5], + [2, 0], + ], + ], + [ + [ + [2, 3], + [3, 3], + [3, 4], + [2, 4], + [2, 3], + ], + ], + ]); + }); + }); +}); diff --git a/packages/cli/src/cog/__test__/cog.vrt.test.ts b/packages/cli/src/cog/__test__/cog.vrt.test.ts index 2d727c8d5e..9a1003ccb0 100644 --- a/packages/cli/src/cog/__test__/cog.vrt.test.ts +++ b/packages/cli/src/cog/__test__/cog.vrt.test.ts @@ -1,13 +1,12 @@ import { EpsgCode } from '@basemaps/geo'; import { FileOperatorSimple, LogConfig, ProjectionTileMatrixSet } from '@basemaps/shared'; +import { qkToName } from '@basemaps/shared/build/tms/__test__/test.util'; import { round } from '@basemaps/test/build/rounding'; -import { FeatureCollection } from 'geojson'; import o from 'ospec'; import { GdalCogBuilder } from '../../gdal/gdal'; -import { Cutline } from '../cutline'; import { CogVrt } from '../cog.vrt'; +import { Cutline } from '../cutline'; import { SourceTiffTestHelper } from './source.tiff.testhelper'; -import { qkToName } from '@basemaps/shared/build/tms/__test__/test.util'; o.spec('cog.vrt', () => { const tmpFolder = '/tmp/my-tmp-folder'; @@ -21,14 +20,9 @@ o.spec('cog.vrt', () => { const googleProj = ProjectionTileMatrixSet.get(EpsgCode.Google); - const [tif1Path, tif2Path] = [1, 2].map((i) => `${testDir}/tif${i}.tiff`); - - const [tif1Poly, tif2Poly] = SourceTiffTestHelper.tiffPolygons(); - - const sourceGeo = { - type: 'FeatureCollection', - features: [], - } as FeatureCollection; + const sourceBounds = SourceTiffTestHelper.tiffNztmBounds(testDir); + const [tif1, tif2] = sourceBounds; + const [tif1Path, tif2Path] = sourceBounds.map(({ name }) => name); let cutTiffArgs: Array> = []; @@ -46,15 +40,12 @@ o.spec('cog.vrt', () => { o.beforeEach(() => { runSpy = o.spy(); - job.projection = job.source.projection = EpsgCode.Google; + job.projection = EpsgCode.Google; + job.source.projection = EpsgCode.Nztm2000; job.source.resZoom = 13; - job.source.files = [tif1Path, tif2Path]; gdal = { run: runSpy }; (GdalCogBuilder as any).getGdal = (): any => gdal; - sourceGeo.features = SourceTiffTestHelper.makeTiffFeatureCollection( - [tif1Poly, tif2Poly], - [tif1Path, tif2Path], - ).features; + job.source.files = [tif1, tif2]; cutTiffArgs = []; FileOperatorSimple.writeJson = ((...args: any): any => { @@ -71,17 +62,17 @@ o.spec('cog.vrt', () => { job.source.resZoom = 17; - const vrt = await CogVrt.buildVrt(tmpFolder, job, sourceGeo, cutline, qkToName('31133322'), logger); + const vrt = await CogVrt.buildVrt(tmpFolder, job, cutline, qkToName('31133322'), logger); - o(job.source.files).deepEquals([tif1Path, tif2Path]); - o(cutline.clipPoly.length).equals(3); + o(job.source.files).deepEquals([tif1, tif2]); + o(cutline.clipPoly.length).equals(2); o(vrt).equals('/tmp/my-tmp-folder/cog.vrt'); }); o('not within tile', async () => { const cutline = new Cutline(googleProj, await Cutline.loadCutline(testDir + '/kapiti.geojson')); - const vrt = await CogVrt.buildVrt(tmpFolder, job, sourceGeo, cutline, qkToName('3131110001'), logger); + const vrt = await CogVrt.buildVrt(tmpFolder, job, cutline, qkToName('3131110001'), logger); o(cutTiffArgs.length).equals(0); o(vrt).equals(null); @@ -89,9 +80,9 @@ o.spec('cog.vrt', () => { }); o('no cutline same projection', async () => { - const vrt = await CogVrt.buildVrt(tmpFolder, job, sourceGeo, new Cutline(googleProj), qkToName('31'), logger); + const vrt = await CogVrt.buildVrt(tmpFolder, job, new Cutline(googleProj), qkToName('31'), logger); - o(job.source.files).deepEquals([tif1Path, tif2Path]); + o(job.source.files).deepEquals([tif1, tif2]); o(cutTiffArgs.length).equals(0); o(vrt).equals('/tmp/my-tmp-folder/cog.vrt'); o(runSpy.callCount).equals(2); @@ -100,9 +91,9 @@ o.spec('cog.vrt', () => { o('no cutline diff projection', async () => { job.projection = EpsgCode.Nztm2000; - const vrt = await CogVrt.buildVrt(tmpFolder, job, sourceGeo, new Cutline(googleProj), qkToName('31'), logger); + const vrt = await CogVrt.buildVrt(tmpFolder, job, new Cutline(googleProj), qkToName('31'), logger); - o(job.source.files).deepEquals([tif1Path, tif2Path]); + o(job.source.files).deepEquals([tif1, tif2]); o(cutTiffArgs.length).equals(0); o(vrt).equals('/tmp/my-tmp-folder/cog.vrt'); o(runSpy.callCount).equals(2); @@ -114,10 +105,10 @@ o.spec('cog.vrt', () => { const name = qkToName('311333222321113310'); - const vrt = await CogVrt.buildVrt(tmpFolder, job, sourceGeo, cutline, name, logger); + const vrt = await CogVrt.buildVrt(tmpFolder, job, cutline, name, logger); o(vrt).equals('/tmp/my-tmp-folder/cog.vrt'); - o(job.source.files).deepEquals([tif2Path]); + o(job.source.files).deepEquals([tif2]); o(cutTiffArgs.length).equals(0); o(runSpy.callCount).equals(2); o(runSpy.args[0]).equals('gdalwarp'); @@ -125,21 +116,50 @@ o.spec('cog.vrt', () => { o('intersected cutline', async () => { const cutline = new Cutline(googleProj, await Cutline.loadCutline(testDir + '/kapiti.geojson'), 20); - - const poly = round(cutline.toGeoJson(), 6); - job.output.cutline = { blend: 20, source: 'cutline.json' }; const name = qkToName('311333222321113'); - const vrt = await CogVrt.buildVrt(tmpFolder, job, sourceGeo, cutline, name, logger); + const vrt = await CogVrt.buildVrt(tmpFolder, job, cutline, name, logger); o(vrt).equals('/tmp/my-tmp-folder/cog.vrt'); - o(job.source.files).deepEquals([tif2Path]); + o(job.source.files).deepEquals([tif2]); o(cutTiffArgs.length).equals(1); o(cutTiffArgs[0][1]).deepEquals(cutline.toGeoJson()); - o(round(cutTiffArgs[0], 6)).deepEquals(['/tmp/my-tmp-folder/cutline.geojson', poly]); + o(round(cutTiffArgs[0], 6)).deepEquals([ + '/tmp/my-tmp-folder/cutline.geojson', + { + type: 'FeatureCollection', + features: [ + { + type: 'Feature', + geometry: { + type: 'MultiPolygon', + coordinates: [ + [ + [ + [174.872663, -40.879256], + [174.902917, -40.879256], + [174.903924, -40.878474], + [174.909682, -40.871987], + [174.919346, -40.866023], + [174.922943, -40.861803], + [174.922943, -40.839788], + [174.913543, -40.839788], + [174.910775, -40.844398], + [174.891215, -40.858532], + [174.879634, -40.870215], + [174.872663, -40.879256], + ], + ], + ], + }, + properties: {}, + }, + ], + }, + ]); }); o('1 surrounded with s3 files', async () => { @@ -151,16 +171,22 @@ o.spec('cog.vrt', () => { const vtif1 = '/vsis3' + tif1Path; const vtif2 = '/vsis3' + tif2Path; - job.source.files = [s3tif1, s3tif2]; + job.source.files = job.source.files.map((o) => { + o = Object.assign({}, o); + o.name = 's3:/' + o.name; + return o; + }); job.output.cutline = { blend: 10, source: 'cutline.json' }; const cutline = new Cutline(googleProj, await Cutline.loadCutline(testDir + '/mana.geojson')); - const vrt = await CogVrt.buildVrt(tmpFolder, job, sourceGeo, cutline, qkToName('3131110001'), logger); + const vrt = await CogVrt.buildVrt(tmpFolder, job, cutline, qkToName('3131110001'), logger); o(vrt).equals('/tmp/my-tmp-folder/cog.vrt'); o(cutTiffArgs.length).equals(1); o(cutTiffArgs[0][0]).equals(tmpFolder + '/cutline.geojson'); + o(cutline.clipPoly.length).equals(1); + const geo = cutline.toGeoJson(); o(geo.type).equals('FeatureCollection'); @@ -179,6 +205,31 @@ o.spec('cog.vrt', () => { ]); } + o(round(coordinates, 5)).deepEquals([ + [ + [ + [174.77125, -41.09044], + [174.77212, -41.09741], + [174.78832, -41.08499], + [174.78135, -41.07763], + [174.78125, -41.07761], + [174.78117, -41.07752], + [174.78108, -41.07748], + [174.78099, -41.07749], + [174.78093, -41.07752], + [174.7809, -41.07756], + [174.78091, -41.07759], + [174.78097, -41.07765], + [174.78097, -41.07771], + [174.78084, -41.07778], + [174.78069, -41.07779], + [174.7801, -41.07774], + [174.77915, -41.08119], + [174.77125, -41.09044], + ], + ], + ]); + o(runSpy.callCount).equals(2); o(mount.calls.map((c: any) => c.args[0])).deepEquals([tmpFolder, s3tif1, s3tif2]); @@ -208,7 +259,7 @@ o.spec('cog.vrt', () => { '-wo', 'NUM_THREADS=ALL_CPUS', '-s_srs', - 'EPSG:3857', + 'EPSG:2193', '-t_srs', 'EPSG:3857', '-tr', @@ -224,7 +275,5 @@ o.spec('cog.vrt', () => { ], logger, ]); - - o(job.source.files).deepEquals([s3tif1, s3tif2]); }); }); diff --git a/packages/cli/src/cog/__test__/cutline.test.ts b/packages/cli/src/cog/__test__/cutline.test.ts index a0eeb8ee42..f76635546d 100644 --- a/packages/cli/src/cog/__test__/cutline.test.ts +++ b/packages/cli/src/cog/__test__/cutline.test.ts @@ -1,4 +1,4 @@ -import { EpsgCode, GeoJson } from '@basemaps/geo'; +import { EpsgCode } from '@basemaps/geo'; import { ProjectionTileMatrixSet } from '@basemaps/shared'; import { Approx } from '@basemaps/test'; import { round } from '@basemaps/test/build/rounding'; @@ -7,10 +7,31 @@ 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'; o.spec('cutline', () => { const testDir = `${__dirname}/../../../__test.assets__`; const googleProj = ProjectionTileMatrixSet.get(EpsgCode.Google); + const nztmProj = 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 name = qkToName('311333222321113310'); + + const job = SourceTiffTestHelper.makeCogJob(); + job.projection = EpsgCode.Google; + job.source.projection = EpsgCode.Nztm2000; + const sourceBounds = SourceTiffTestHelper.tiffNztmBounds(testDir); + const [tif1, tif2] = sourceBounds; + job.source.files = [tif1, tif2]; + + cutline.filterSourcesForName(name, job); + + o(job.source.files).deepEquals([tif2]); + }); + }); o('loadCutline', async () => { const cutline = new Cutline(googleProj, await Cutline.loadCutline(testDir + '/mana.geojson')); @@ -38,24 +59,26 @@ o.spec('cutline', () => { o('findCovering', async () => { const cutline = new Cutline(googleProj, await Cutline.loadCutline(testDir + '/mana.geojson')); - const feature = SourceTiffTestHelper.makeTiffFeatureCollection(); + const bounds = SourceTiffTestHelper.tiffNztmBounds(); o(cutline.clipPoly.length).equals(2); - const result = (cutline as any).findCovering({ bounds: feature, resZoom: 15 }); + (cutline as any).findCovering({ projection: EpsgCode.Nztm2000, bounds, resZoom: 15 }); - o(result.length).equals(1); + o((cutline as any).srcPoly.length).equals(1); }); o.spec('optmize', async () => { - const bounds = SourceTiffTestHelper.makeTiffFeatureCollection(); + const bounds = SourceTiffTestHelper.tiffNztmBounds(); o('nztm', () => { - const proj = ProjectionTileMatrixSet.get(EpsgCode.Nztm2000); + const cutline = new Cutline(nztmProj); - const cutline = new Cutline(proj); - - const covering = cutline.optimizeCovering({ bounds, resZoom: 14 } as SourceMetadata); + const covering = cutline.optimizeCovering({ + projection: EpsgCode.Nztm2000, + bounds, + resZoom: 14, + } as SourceMetadata); o(covering.length).equals(covering.length); o(covering[4]).deepEquals({ @@ -81,23 +104,86 @@ o.spec('cutline', () => { ]); }); - o('boundary part inland, part coastal', async () => { - const poly = GeoJson.toPositionPolygon([174.89, -40.83, 174.92, -40.85]); + o('source polygon is smoothed', async () => { + const tiff1 = { + name: '/path/to/tiff1.tiff', + x: 1759315.9336568, + y: 5476120.089572778, + width: 2577.636033663526, + height: 2275.3233966259286, + }; + + const tiff2 = { + name: '/path/to/tiff2.tiff', + x: 1759268.0010635862, + y: 5473899.771969615, + width: 2576.894309356343, + height: 2275.3359155077487, + }; + + const bounds = [tiff1, tiff2]; + const cutline = new Cutline(googleProj, await Cutline.loadCutline(testDir + '/kapiti.geojson'), 500); - const bounds = SourceTiffTestHelper.makeTiffFeatureCollection(poly); + const covering = cutline.optimizeCovering({ + projection: EpsgCode.Nztm2000, + bounds, + resZoom: 22, + } as SourceMetadata); + + o(round(cutline.clipPoly, 4)).deepEquals([ + [ + [ + [19468580.0838, -4993280.8853], + [19471098.3762, -4993280.8853], + [19471932.5568, -4992600.0637], + [19472092.0139, -4992352.6943], + [19472092.0139, -4987203.6846], + [19471983.8261, -4987203.6846], + [19470978.4379, -4989417.4454], + [19468800.9775, -4991497.5405], + [19468580.0838, -4991792.2036], + [19468580.0838, -4993280.8853], + ], + ], + ]); + + o(covering.map((c) => c.name)).deepEquals([ + '14-16151-10232', + '14-16151-10233', + '14-16152-10231', + '14-16152-10232', + '14-16152-10233', + '18-258444-163695', + ]); + }); + + o('boundary part inland, part coastal', async () => { const cutline = new Cutline(googleProj, await Cutline.loadCutline(testDir + '/kapiti.geojson'), 500); + const bounds = [ + { + name: 'tiff1', + x: 1759315.9336568, + y: 5476120.089572778, + width: 2577.636033663526, + height: 2275.3233966259286, + }, + ]; - const covering = cutline.optimizeCovering({ bounds, resZoom: 22 } as SourceMetadata); + const covering = cutline.optimizeCovering({ + projection: EpsgCode.Nztm2000, + bounds, + resZoom: 22, + } as SourceMetadata); o(round(cutline.clipPoly, 4)).deepEquals([ [ [ - [19470094.8286, -4990261.5441], - [19472027.7232, -4990261.5441], - [19472027.7232, -4987279.214], - [19471949.5241, -4987279.214], + [19470015.7313, -4990337.1045], + [19472091.9684, -4990337.1045], + [19472091.9684, -4987203.6846], + [19471983.8261, -4987203.6846], [19470978.4379, -4989417.4454], - [19470094.8286, -4990261.5441], + [19470015.7313, -4990337.1045], ], ], ]); @@ -106,14 +192,18 @@ o.spec('cutline', () => { '14-16152-10231', '15-32304-20464', '15-32305-20464', - '17-129222-81847', + '18-258444-163695', ]); }); o('low res', () => { const cutline = new Cutline(googleProj); - const covering = cutline.optimizeCovering({ bounds, resZoom: 13 } as SourceMetadata); + const covering = cutline.optimizeCovering({ + projection: EpsgCode.Nztm2000, + bounds, + resZoom: 13, + } as SourceMetadata); o(covering.length).equals(covering.length); o(covering.map((c) => c.name)).deepEquals(['8-252-159', '8-252-160']); @@ -121,6 +211,7 @@ o.spec('cutline', () => { o('hi res', () => { const covering2 = new Cutline(googleProj).optimizeCovering({ + projection: EpsgCode.Nztm2000, bounds, resZoom: 18, } as SourceMetadata); diff --git a/packages/cli/src/cog/__test__/source.tiff.testhelper.ts b/packages/cli/src/cog/__test__/source.tiff.testhelper.ts index 9dbcd0823f..3e61436d52 100644 --- a/packages/cli/src/cog/__test__/source.tiff.testhelper.ts +++ b/packages/cli/src/cog/__test__/source.tiff.testhelper.ts @@ -1,4 +1,5 @@ -import { GeoJson, EpsgCode } from '@basemaps/geo'; +import { EpsgCode } from '@basemaps/geo'; +import { NamedBounds } from '@basemaps/shared'; import { CogJob } from '../types'; export const SourceTiffTestHelper = { @@ -6,7 +7,7 @@ export const SourceTiffTestHelper = { return { projection: EpsgCode.Google, source: { - files: [] as string[], + files: [] as NamedBounds[], resZoom: 13, pixelScale: 9.55, path: '', @@ -20,35 +21,26 @@ export const SourceTiffTestHelper = { } as CogJob; }, - tiffPolygons(): GeoJSON.Position[][] { + tiffNztmBounds(path = '/path/to'): NamedBounds[] { return [ - [ - [174.56568549279925, -40.83842049282911], - [174.57337757492053, -41.162595138463644], - [174.85931317513828, -41.15833331417625], - [174.85022498195937, -40.834206771481526], - [174.56568549279925, -40.83842049282911], - ], - - [ - [174.85022498195937, -40.834206771481526], - [174.85931317513828, -41.15833331417625], - [175.14518230301843, -41.15336229204068], - [175.13469922175867, -40.829291849263555], - [174.85022498195937, -40.834206771481526], - ], + { + name: path + '/tiff1.tiff', + x: 1732000, + y: 5442000, + width: 24000, + height: 36000, + }, + { + name: path + '/tiff2.tiff', + x: 1756000, + y: 5442000, + width: 24000, + height: 36000, + }, ]; }, - polygonsToPaths(polygons: GeoJSON.Position[][]): string[] { - return polygons.map((_, i): string => `/path/to/tiff${i}.tiff`); - }, - - makeTiffFeatureCollection(polygons?: GeoJSON.Position[][], paths?: string[]): GeoJSON.FeatureCollection { - if (polygons == null) polygons = SourceTiffTestHelper.tiffPolygons(); - if (paths == null) paths = SourceTiffTestHelper.polygonsToPaths(polygons); - return GeoJson.toFeatureCollection( - polygons.map((coords, i) => GeoJson.toFeaturePolygon([coords], { tiff: paths![i] })), - ); + namedBoundsToPaths(bounds: NamedBounds[]): string[] { + return bounds.map(({ name }): string => `/path/to/tiff${name}.tiff`); }, }; diff --git a/packages/cli/src/cog/builder.ts b/packages/cli/src/cog/builder.ts index a473c8d27a..449c175822 100644 --- a/packages/cli/src/cog/builder.ts +++ b/packages/cli/src/cog/builder.ts @@ -1,4 +1,4 @@ -import { Epsg, GeoJson, Bounds } from '@basemaps/geo'; +import { Bounds, Epsg } from '@basemaps/geo'; import { CompositeError, FileOperatorSimple, @@ -16,7 +16,7 @@ import { Cutline } from './cutline'; // import { CogBuilderMetadata, SourceMetadata } from './types'; export const InvalidProjectionCode = 32767; -export const CacheVersion = 'v2'; // bump this number to invalidate the cache +export const CacheVersion = 'v3'; // bump this number to invalidate the cache export const CacheFolder = './.cache'; /** @@ -69,7 +69,8 @@ export class CogBuilder { let projection = this.srcProj; let nodata: number | undefined; let count = 0; - const coordinates = sources.map((source) => { + + const queue = sources.map((source) => { return this.q(async () => { count++; if (count % 50 == 0) { @@ -86,7 +87,10 @@ export class CogBuilder { bands = tiffBandCount.length; } - const output = this.getTifBounds(tiff); + const name = source instanceof CogSourceFile ? source.fileName : source.name; + + const output = { ...Bounds.fromBbox(image.bbox).toJson(), name }; + if (CogSourceFile.isSource(source)) { await source.close(); } @@ -120,7 +124,7 @@ export class CogBuilder { }); }); - const polygons = await Promise.all(coordinates); + const bounds = await Promise.all(queue); if (projection == null) throw new Error('No projection detected'); if (resX == -1) throw new Error('No resolution detected'); @@ -130,7 +134,7 @@ export class CogBuilder { projection: projection.code, nodata, bands, - bounds: GeoJson.toFeatureCollection(polygons), + bounds, pixelScale: resX, resZoom: this.targetProj.getTiffResZoom(resX), }; @@ -185,34 +189,6 @@ export class CogBuilder { return noDataNum; } - /** - * Generate the bounding boxes for a GeoTiff converting to WGS84 - * @param tiff - */ - getTifBounds(tiff: CogTiff): GeoJSON.Feature { - const image = tiff.getImage(0); - const bbox = image.bbox; - const topLeft = [bbox[0], bbox[3]]; - const topRight = [bbox[2], bbox[3]]; - const bottomRight = [bbox[2], bbox[1]]; - const bottomLeft = [bbox[0], bbox[1]]; - - const projection = this.findProjection(tiff); - const projProjection = ProjectionTileMatrixSet.tryGet(projection.code); - if (projProjection == null) { - this.logger.error({ tiff: tiff.source.name }, 'Failed to get tiff projection'); - throw new Error('Invalid tiff projection: ' + projection); - } - - const { toWsg84 } = projProjection; - - const points = [ - [toWsg84(topLeft), toWsg84(bottomLeft), toWsg84(bottomRight), toWsg84(topRight), toWsg84(topLeft)], - ]; - - return GeoJson.toFeaturePolygon(points, { tiff: tiff.source.name }); - } - /** Cache the bounds lookup so we do not have to requery the bounds between CLI calls */ private async getMetadata(tiffs: CogSource[]): Promise { const cacheKey = diff --git a/packages/cli/src/cog/clipped.multipolygon.ts b/packages/cli/src/cog/clipped.multipolygon.ts new file mode 100644 index 0000000000..12cafadf6e --- /dev/null +++ b/packages/cli/src/cog/clipped.multipolygon.ts @@ -0,0 +1,31 @@ +import { Bounds } from '@basemaps/geo'; +import lineclip from 'lineclip'; +import pc, { MultiPolygon } from 'polygon-clipping'; + +function samePoint(a: number[], b: number[]): boolean { + return a[0] == b[0] && a[1] == b[1]; +} + +export function clipMultipolygon(polygons: MultiPolygon, bounds: Bounds): MultiPolygon { + const result: MultiPolygon = []; + const bbox = bounds.toBbox(); + for (const poly of polygons) { + const clipped = lineclip.polygon(poly[0], bbox); + if (clipped.length != 0) { + if (!samePoint(clipped[0], clipped[clipped.length - 1])) clipped.push(clipped[0]); + result.push([clipped]); + } + } + return result; +} + +export function removeDegenerateEdges(polygons: MultiPolygon, bounds: Bounds): MultiPolygon { + return pc.intersection(polygons, bounds.toPolygon()); +} + +export function polyContainsBounds(poly: MultiPolygon, bounds: Bounds): boolean { + const clipped = clipMultipolygon(poly, bounds); + if (clipped.length != 1 || clipped[0].length != 1 || clipped[0][0].length != 5) return false; + + return Bounds.fromMultiPolygon(clipped).containsBounds(bounds); +} diff --git a/packages/cli/src/cog/cog.vrt.ts b/packages/cli/src/cog/cog.vrt.ts index 5457fd2b15..9cebfffeba 100644 --- a/packages/cli/src/cog/cog.vrt.ts +++ b/packages/cli/src/cog/cog.vrt.ts @@ -1,11 +1,10 @@ -import { FileOperator, LogType, isConfigS3Role, Aws } from '@basemaps/shared'; -import { FeatureCollection } from 'geojson'; -import { Cutline } from './cutline'; -import { CogJob } from './types'; -import { onProgress } from './cog'; -import { GdalCogBuilder } from '../gdal/gdal'; import { Epsg } from '@basemaps/geo'; +import { Aws, FileOperator, isConfigS3Role, LogType } from '@basemaps/shared'; +import { GdalCogBuilder } from '../gdal/gdal'; import { GdalCommand } from '../gdal/gdal.command'; +import { onProgress } from './cog'; +import { Cutline } from './cutline'; +import { CogJob } from './types'; /** * Build the VRT for the needed source imagery @@ -28,14 +27,14 @@ async function buildPlainVrt( gdalCommand.setCredentials(credentials); } else { if (gdalCommand.mount) { - for (const file of job.source.files) { - gdalCommand.mount(file); + for (const { name } of job.source.files) { + gdalCommand.mount(name); } } } logger.debug({ buildOpts: buildOpts.join(' ') }, 'gdalbuildvrt'); - const sourceFiles = job.source.files.map((c) => c.replace('s3://', '/vsis3/')); + const sourceFiles = job.source.files.map(({ name }) => name.replace('s3://', '/vsis3/')); await gdalCommand.run('gdalbuildvrt', [...buildOpts, vrtPath, ...sourceFiles], logger); } @@ -79,13 +78,12 @@ async function buildWarpVrt( export const CogVrt = { /** - * Build a vrt file for a `name` set with some tiffs transformed with a cutline + * Build a vrt file for a COG `name` that transforms the source imagery with a cutline * * @param tmpFolder temporary `vrt` and `cutline.geojson` will be written here * @param job - * @param sourceGeo a GeoJSON object which contains the boundaries for the source imagery - * @param cutline Used to filter the sources and cutline - * @param name to reduce vrt and cutline + * @param cutline Used to filter the source imagery + * @param name COG tile to reduce vrt and cutline * @param logger * * @return the path to the vrt file @@ -93,7 +91,6 @@ export const CogVrt = { async buildVrt( tmpFolder: string, job: CogJob, - sourceGeo: FeatureCollection, cutline: Cutline, name: string, logger: LogType, @@ -102,7 +99,7 @@ export const CogVrt = { const inputTotal = job.source.files.length; - cutline.filterSourcesForName(name, job, sourceGeo); + cutline.filterSourcesForName(name, job); if (job.source.files.length == 0) { return null; diff --git a/packages/cli/src/cog/cutline.ts b/packages/cli/src/cog/cutline.ts index e3399fb75c..452c9c95d0 100644 --- a/packages/cli/src/cog/cutline.ts +++ b/packages/cli/src/cog/cutline.ts @@ -1,10 +1,10 @@ import { Bounds, Epsg, GeoJson, Tile, TileMatrixSet } from '@basemaps/geo'; import { compareName, FileOperator, NamedBounds, ProjectionTileMatrixSet } from '@basemaps/shared'; -import { FeatureCollection, Position } from 'geojson'; -import { basename } from 'path'; -import pc, { MultiPolygon, Polygon, Ring } from 'polygon-clipping'; +import { FeatureCollection } from 'geojson'; 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'; const { intersection, union } = pc; /** Padding to always apply to image boundies */ @@ -14,18 +14,13 @@ function findGeoJsonProjection(geojson: any | null): Epsg { return Epsg.parse(geojson?.crs?.properties?.name ?? '') ?? Epsg.Wgs84; } -function polysSameShape(a: Polygon, b: Polygon): boolean { - if (a.length != b.length) return false; - for (let i = 0; i < a.length; ++i) { - if (a[i].length != b[i].length) return false; - } - return true; -} - function namedBounds(tms: TileMatrixSet, tile: Tile): NamedBounds { return { name: TileMatrixSet.tileToName(tile), ...tms.tileToSourceBounds(tile).toJson() }; } +/** + * Filter out duplicate tiles + */ function addNonDupes(list: Tile[], addList: Tile[]): void { const len = list.length; for (const add of addList) { @@ -47,6 +42,7 @@ export class Cutline { targetProj: ProjectionTileMatrixSet; blend: number; tms: TileMatrixSet; // convience to targetProj.tms + private srcPoly: MultiPolygon = []; /** * Create a Cutline instance from a `GeoJSON FeatureCollection`. @@ -98,38 +94,36 @@ export class Cutline { * @param job * @param sourceGeo */ - filterSourcesForName(name: string, job: CogJob, sourceGeo: FeatureCollection): void { + filterSourcesForName(name: string, job: CogJob): void { const tile = TileMatrixSet.nameToTile(name); - const qkBounds = this.targetProj.tileToSourceBounds(tile); + const { targetProj } = this; + const qkBounds = targetProj.tileToSourceBounds(tile); const qkPadded = this.padBounds(qkBounds, job.source.resZoom); - const srcTiffs = new Set(); + let qkBoundsInSrcProj = qkPadded; - for (const f of sourceGeo.features) { - const { geometry } = f; - if (geometry.type === 'Polygon') { - const srcBounds = this.sourceWsg84PolyToBounds(geometry.coordinates[0]); - if (qkPadded.intersects(srcBounds)) { - srcTiffs.add(basename(f.properties!.tiff!)); - } - } + if (job.source.projection !== this.tms.projection.code) { + // convert the padded quadKey to source projection ensuring fully enclosed + const sourceProj = ProjectionTileMatrixSet.get(job.source.projection); + const poly = [qkBoundsInSrcProj.toPolygon()]; + targetProj.projectMultipolygon(poly, sourceProj); + qkBoundsInSrcProj = Bounds.fromMultiPolygon(poly); } - job.source.files = job.source.files.filter((path) => srcTiffs.has(basename(path))); + + job.source.files = job.source.files.filter((image) => qkBoundsInSrcProj.intersects(Bounds.fromJson(image))); if (this.clipPoly.length > 0) { - const boundsPoly = GeoJson.toPositionPolygon(qkPadded.toBbox()) as Polygon; - const poly = intersection(boundsPoly, this.clipPoly); - if (poly == null || poly.length == 0) { + const poly = clipMultipolygon(this.clipPoly, qkPadded); + if (poly.length == 0) { // this tile is not needed this.clipPoly = []; job.source.files = []; - } else if ( - poly.length == 1 && - polysSameShape(poly[0], boundsPoly) && - this.sourcePolyToBounds(poly[0][0]).containsBounds(qkBounds) - ) { + } else if (polyContainsBounds(poly, qkBounds)) { // tile is completely surrounded; no cutline polygons needed this.clipPoly = []; + } else { + // set the cutline polygons to just the area of interest (minus degenerate edges) + this.clipPoly = removeDegenerateEdges(poly, qkPadded); } } } @@ -139,7 +133,7 @@ export class Cutline { * @param featureCollection Source TIff Polygons in GeoJson WGS84 */ optimizeCovering(sourceMetadata: SourceMetadata): NamedBounds[] { - const srcArea = this.findCovering(sourceMetadata); + this.findCovering(sourceMetadata); const { resZoom } = sourceMetadata; // Look for the biggest tile size we are allowed to create. @@ -154,7 +148,7 @@ export class Cutline { for (const tile of tms.topLevelTiles()) { // Don't make COGs with a tile.z < minZ. - tiles = tiles.concat(this.makeTiles(tile, srcArea, minZ, CoveringFraction).tiles); + tiles = tiles.concat(this.makeTiles(tile, this.srcPoly, minZ, CoveringFraction).tiles); } if (tiles.length == 0) { @@ -181,36 +175,6 @@ export class Cutline { return GeoJson.toFeatureCollection([GeoJson.toFeatureMultiPolygon(clipPoly.map((p) => [p[0].map(toWsg84)]))]); } - /** - * Convert a Wsg84 "square" polygon to it's bounds - */ - private sourcePolyToBounds(poly: Position[]): Bounds { - const [x1, y1] = poly[0]; - const [x2, y2] = poly[2]; - - return new Bounds(Math.min(x1, x2), Math.min(y1, y2), Math.abs(x2 - x1), Math.abs(y2 - y1)); - } - - /** - * Convert a Wsg84 polygon to it's bounds - */ - private sourceWsg84PolyToBounds(poly: Position[]): Bounds { - const { fromWsg84 } = this.targetProj; - let minX = Infinity; - let minY = Infinity; - let maxX = -Infinity; - let maxY = -Infinity; - for (const p of poly) { - const [x, y] = fromWsg84(p); - if (x < minX) minX = x; - if (y < minY) minY = y; - if (x > maxX) maxX = x; - if (y > maxY) maxY = y; - } - - return new Bounds(minX, minY, maxX - minX, maxY - minY); - } - /** * Merge child nodes that have at least a covering fraction * @param tile the tile to descend @@ -226,12 +190,14 @@ export class Cutline { minZ: number, coveringFraction: number, ): { tiles: Tile[]; fractionCovered: number } { - const clip = this.targetProj.tileToPolygon(tile) as Polygon; - const intArea = intersection(srcArea, clip); + const clipBounds = this.targetProj.tileToSourceBounds(tile); - if (intArea.length == 0) { + srcArea = clipMultipolygon(srcArea, clipBounds); + + if (srcArea.length == 0) { return { tiles: [], fractionCovered: 0 }; } + if (tile.z == minZ + 4) { return { tiles: [tile], fractionCovered: 1 }; } @@ -239,7 +205,7 @@ export class Cutline { const ans = { tiles: [] as Tile[], fractionCovered: 0 }; for (const child of this.tms.coverTile(tile)) { - const { tiles, fractionCovered } = this.makeTiles(child, intArea, minZ, coveringFraction); + const { tiles, fractionCovered } = this.makeTiles(child, srcArea, minZ, coveringFraction); if (fractionCovered != 0) { ans.fractionCovered += fractionCovered * 0.25; addNonDupes(ans.tiles, tiles); @@ -266,33 +232,38 @@ export class Cutline { * @param sourceMetadata */ - private findCovering(sourceMetadata: SourceMetadata): MultiPolygon { - let srcArea: MultiPolygon = []; + private findCovering(sourceMetadata: SourceMetadata): void { + let srcPoly: MultiPolygon = []; const { resZoom } = sourceMetadata; // merge imagery bounds - for (const { geometry } of sourceMetadata.bounds.features) { - if (geometry.type === 'Polygon') { - // ensure source polys overlap by using their bounding box - const poly = GeoJson.toPositionPolygon( - this.padBounds(this.sourceWsg84PolyToBounds(geometry.coordinates[0]), resZoom).toBbox(), - ) as Polygon; - if (srcArea.length == 0) { - srcArea.push(poly); - } else { - srcArea = union(srcArea, poly) as MultiPolygon; - } - } + for (const image of sourceMetadata.bounds) { + const poly = Bounds.fromJson(image).toPolygon() as Polygon; + srcPoly = union(srcPoly, poly); } - if (this.clipPoly.length != 0) { - // clip the imagery bounds - this.clipPoly = (intersection(srcArea, this.clipPoly) ?? []) as MultiPolygon; - return this.clipPoly; - } + // Convert polygon to target projection + const sourceProj = ProjectionTileMatrixSet.get(sourceMetadata.projection); + sourceProj.projectMultipolygon(srcPoly, this.targetProj); + this.srcPoly = srcPoly; + + if (this.clipPoly.length == 0) return; - // no cutline return imagery bounds - return srcArea; + const srcBounds = Bounds.fromMultiPolygon(srcPoly); + const boundsPadded = this.padBounds(srcBounds, resZoom); + + const poly = clipMultipolygon(this.clipPoly, boundsPadded); + if (poly == null || poly.length == 0) { + throw new Error('No intersection between source imagery and cutline'); + } + if (polyContainsBounds(poly, srcBounds)) { + // tile is completely surrounded; no cutline polygons needed + this.clipPoly = []; + } else { + // set the cutline polygons to just the area of interest (minus degenerate edges) + this.clipPoly = removeDegenerateEdges(poly, boundsPadded); + this.srcPoly = intersection(srcPoly, this.clipPoly) ?? []; + } } /** @@ -305,6 +276,8 @@ export class Cutline { const px = this.tms.pixelScale(resZoom); // Ensure cutline blend does not interferre with non-costal edges - return bounds.scaleFromCenter((bounds.width + px * (PixelPadding + this.blend) * 2) / bounds.width); + const widthScale = (bounds.width + px * (PixelPadding + this.blend) * 2) / bounds.width; + const heightScale = (bounds.height + px * (PixelPadding + this.blend) * 2) / bounds.height; + return bounds.scaleFromCenter(widthScale, heightScale); } } diff --git a/packages/cli/src/cog/job.ts b/packages/cli/src/cog/job.ts index d2f06cad35..b47035e875 100644 --- a/packages/cli/src/cog/job.ts +++ b/packages/cli/src/cog/job.ts @@ -123,6 +123,8 @@ 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 { tms } = cutline.targetProj; const files = metadata.files.sort(Bounds.compareArea); @@ -180,7 +182,7 @@ export const CogJobFactory = { pixelScale: metadata.pixelScale, resZoom: metadata.resZoom, projection: metadata.projection, - files: tiffList, + files: metadata.bounds, options: { maxConcurrency }, }, bounds: metadata.targetBounds, @@ -201,13 +203,13 @@ export const CogJobFactory = { const jobFile = getJobPath(job, `job.json`); await outputFs.writeJson(jobFile, job, logger); - if (ctx.cutline != null) { + if (cutline.clipPoly.length != 0) { const geoJsonCutlineOutput = getJobPath(job, `cutline.geojson.gz`); await outputFs.writeJson(geoJsonCutlineOutput, cutline.toGeoJson(), logger); } const geoJsonSourceOutput = getJobPath(job, `source.geojson`); - await outputFs.writeJson(geoJsonSourceOutput, metadata.bounds, logger); + 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); diff --git a/packages/cli/src/cog/types.ts b/packages/cli/src/cog/types.ts index 148a87d977..5cb5b2d57a 100644 --- a/packages/cli/src/cog/types.ts +++ b/packages/cli/src/cog/types.ts @@ -14,7 +14,7 @@ export interface CogJob { source: { /** List of input files */ - files: string[]; + files: NamedBounds[]; /** The number of pixels per meter for the best source image */ pixelScale: number; @@ -78,7 +78,7 @@ export interface SourceMetadata { /** Number of imagery bands generally RGB (3) or RGBA (4) */ bands: number; /** Bounding box for polygons */ - bounds: GeoJSON.FeatureCollection; + bounds: NamedBounds[]; /** The number of pixels per meter for the best source image */ pixelScale: number; diff --git a/packages/cli/src/lineclip.d.ts b/packages/cli/src/lineclip.d.ts new file mode 100644 index 0000000000..f67ebc3d91 --- /dev/null +++ b/packages/cli/src/lineclip.d.ts @@ -0,0 +1,3 @@ +declare module 'lineclip' { + export function polygon(points: number[][], bbox: number[]): [number, number][]; +} diff --git a/packages/geo/src/__tests__/bounds.test.ts b/packages/geo/src/__tests__/bounds.test.ts index e9c4b66de3..c6b2036ab8 100644 --- a/packages/geo/src/__tests__/bounds.test.ts +++ b/packages/geo/src/__tests__/bounds.test.ts @@ -1,6 +1,6 @@ +import { Approx } from '@basemaps/test'; import o from 'ospec'; import { Bounds } from '../bounds'; -import { Approx } from '@basemaps/test'; const TILE_SIZE = 256; function getTile(x = 0, y = 0): Bounds { @@ -110,4 +110,32 @@ o.spec('Bounds', () => { o(compareArea(new Bounds(3, 5, 1, 1), new Bounds(3, 5, 2, 5))).equals(-9); }); + + o('fromMultiPolygon', () => { + const poly = [ + [ + [ + [84, -49], + [74, -40], + [90, -57], + [94, -39], + [84, -49], + ], + ], + ]; + + o(Bounds.fromMultiPolygon(poly).toJson()).deepEquals({ x: 74, y: -57, width: 20, height: 18 }); + }); + + o('toPolygon', () => { + o(new Bounds(74, -57, 20, 18).toPolygon()).deepEquals([ + [ + [74, -57], + [74, -39], + [94, -39], + [94, -57], + [74, -57], + ], + ]); + }); }); diff --git a/packages/geo/src/bounds.ts b/packages/geo/src/bounds.ts index e558c2d928..d7af925c83 100644 --- a/packages/geo/src/bounds.ts +++ b/packages/geo/src/bounds.ts @@ -114,6 +114,21 @@ export class Bounds implements BoundingBox { return [this.x, this.y, this.right, this.bottom]; } + /** + * Convert to GeoJson Polygon coordinates + */ + public toPolygon(): [number, number][][] { + return [ + [ + [this.x, this.y], + [this.x, this.bottom], + [this.right, this.bottom], + [this.right, this.y], + [this.x, this.y], + ], + ]; + } + /** * Scale the bounding box adjusting top, left, width & height * @@ -153,6 +168,30 @@ export class Bounds implements BoundingBox { return new Bounds(this.x - bounds.x, this.y - bounds.y, this.width, this.height); } + /** + * Find the bounds of a GeoJson MultiPolygon + + * @param multipoly the polygon to measure + * @throws Error if `multipoly` is empty + */ + public static fromMultiPolygon(multipoly: number[][][][]): Bounds { + if (multipoly.length == 0) return new Bounds(0, 0, 0, 0); + let minX = multipoly[0][0][0][0]; + let minY = multipoly[0][0][0][1]; + let maxX = minX; + let maxY = minY; + for (const [poly] of multipoly) { + for (const [x, y] of poly) { + if (x < minX) minX = x; + else if (x > maxX) maxX = x; + if (y < minY) minY = y; + else if (y > maxY) maxY = y; + } + } + + return new Bounds(minX, minY, maxX - minX, maxY - minY); + } + /** * Convert a BBox(minX, minY, maxX, maxY) to Bounds(x,y, width, height). * Takes into account the antimeridian. diff --git a/packages/geo/src/geo.json.ts b/packages/geo/src/geo.json.ts index 5a16a74d98..258273866f 100644 --- a/packages/geo/src/geo.json.ts +++ b/packages/geo/src/geo.json.ts @@ -12,18 +12,6 @@ export const GeoJson = { }; }, - toPositionPolygon(bbox: [number, number, number, number]): Position[][] { - return [ - [ - [bbox[0], bbox[1]], - [bbox[0], bbox[3]], - [bbox[2], bbox[3]], - [bbox[2], bbox[1]], - [bbox[0], bbox[1]], - ], - ]; - }, - toPolygon(coordinates: Position[][]): Polygon { return { type: 'Polygon', diff --git a/packages/shared/src/tms/__test__/projection.tile.matrix.set.test.ts b/packages/shared/src/tms/__test__/projection.tile.matrix.set.test.ts index 2c8c4fce70..bb85b938de 100644 --- a/packages/shared/src/tms/__test__/projection.tile.matrix.set.test.ts +++ b/packages/shared/src/tms/__test__/projection.tile.matrix.set.test.ts @@ -107,6 +107,31 @@ o.spec('ProjectionTileMatrixSet', () => { ]); }); + o('projectMultipolygon', () => { + const poly = [ + Bounds.fromBbox([ + 18494091.86765497, + -6051366.655280836, + 19986142.659781612, + -4016307.214216303, + ]).toPolygon(), + ]; + + googleProj.projectMultipolygon(poly, nztmProj); + + o(round(poly, 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/projection.tile.matrix.set.ts b/packages/shared/src/tms/projection.tile.matrix.set.ts index 1d52ba1c6c..30b9f7b16d 100644 --- a/packages/shared/src/tms/projection.tile.matrix.set.ts +++ b/packages/shared/src/tms/projection.tile.matrix.set.ts @@ -64,6 +64,20 @@ export class ProjectionTileMatrixSet { return (epsgCode && CodeMap.get(epsgCode)) ?? null; } + /** + * Project the points in a MultiPolygon array to the `targetProjection`. Modifies `multipoly`. + * Does not support holes in polygons. + */ + projectMultipolygon(multipoly: number[][][][], targetProjection: ProjectionTileMatrixSet): void { + if (targetProjection.tms.projection.code === this.tms.projection.code) return; + + const { toWsg84 } = this; + const { fromWsg84 } = targetProjection; + for (const poly of multipoly) { + poly[0] = poly[0].map((p) => fromWsg84(toWsg84(p))); + } + } + /** * Convert source `[x, y]` coordinates to `[lon, lat]` */ @@ -113,31 +127,17 @@ export class ProjectionTileMatrixSet { return { lat, lon }; } - /** - * Convert a tile to a GeoJson Polygon in Source units - */ - tileToPolygon(tile: Tile): Position[][] { - return [ - [ - [tile.x, tile.y + 1], - [tile.x, tile.y], - [tile.x + 1, tile.y], - [tile.x + 1, tile.y + 1], - [tile.x, tile.y + 1], - ].map(([x, y]) => { - const p = this.tms.tileToSource({ x, y, z: tile.z }); - return [p.x, p.y]; - }), - ]; - } - /** 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( - [GeoJson.toPositionPolygon(Bounds.fromJson(bounds).toBbox())[0].map((p) => toWsg84(p))], + [ + Bounds.fromJson(bounds) + .toPolygon()[0] + .map((p) => toWsg84(p)), + ], { name: bounds.name, }, diff --git a/yarn.lock b/yarn.lock index 530ab61471..27eba3dad6 100644 --- a/yarn.lock +++ b/yarn.lock @@ -4658,6 +4658,11 @@ levn@^0.3.0, levn@~0.3.0: prelude-ls "~1.1.2" type-check "~0.3.2" +lineclip@^1.1.5: + version "1.1.5" + resolved "https://registry.yarnpkg.com/lineclip/-/lineclip-1.1.5.tgz#2bf26067d94354feabf91e42768236db5616fd13" + integrity sha1-K/JgZ9lDVP6r+R5CdoI221YW/RM= + load-json-file@^1.0.0: version "1.1.0" resolved "https://registry.yarnpkg.com/load-json-file/-/load-json-file-1.1.0.tgz#956905708d58b4bab4c2261b04f59f31c99374c0"