Skip to content

Commit

Permalink
fix(cli): mitigate polygon intersection errors
Browse files Browse the repository at this point in the history 
Floating-point rounding can cause the cutline polygon to self cross so we look for points too close together.

The Martinez polygon-clipping libraries currently suffer from errors when lines align; try to
mitigate around this by using simple bounds intersection where appropriate. Other intersection
libraries are intractable for the size of coastline polygons.

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.
  • Loading branch information
Geoff Jacobsen committed Jun 28, 2020
1 parent 42ded3a commit 8d2c937
Show file tree
Hide file tree
Showing 15 changed files with 566 additions and 322 deletions.
4 changes: 1 addition & 3 deletions packages/cli/src/cli/cogify/action.cog.ts
View file
Original file line number Diff line number Diff line change
Expand Up @@ -99,8 +99,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');
Expand All @@ -115,7 +113,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');
Expand Down
99 changes: 68 additions & 31 deletions packages/cli/src/cog/__test__/builder.test.ts
View file
Original file line number Diff line number Diff line change
@@ -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';

Expand All @@ -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,
});
});
});
Expand Down
59 changes: 26 additions & 33 deletions packages/cli/src/cog/__test__/cog.vrt.test.ts
View file
Original file line number Diff line number Diff line change
@@ -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';
Expand All @@ -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<Array<any>> = [];

Expand All @@ -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 => {
Expand All @@ -71,27 +62,27 @@ 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(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);
o(runSpy.callCount).equals(0);
});

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);
Expand All @@ -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);
Expand All @@ -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');
Expand All @@ -129,10 +120,10 @@ o.spec('cog.vrt', () => {

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());

Expand Down Expand Up @@ -180,11 +171,15 @@ 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);
Expand Down Expand Up @@ -264,7 +259,7 @@ o.spec('cog.vrt', () => {
'-wo',
'NUM_THREADS=ALL_CPUS',
'-s_srs',
'EPSG:3857',
'EPSG:2193',
'-t_srs',
'EPSG:3857',
'-tr',
Expand All @@ -280,7 +275,5 @@ o.spec('cog.vrt', () => {
],
logger,
]);

o(job.source.files).deepEquals([s3tif1, s3tif2]);
});
});
Loading

0 comments on commit 8d2c937

Please sign in to comment.