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 03f8f046fc..e4cf21c429 100644 --- a/yarn.lock +++ b/yarn.lock @@ -4613,6 +4613,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"