Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

fix(shared): handle bounds crossing antimeridian #925

Merged
merged 2 commits into from
Jul 16, 2020
Merged

fix(shared): handle bounds crossing antimeridian #925

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6c83260
fix(shared): handle bounds crossing antimeridian
Jul 14, 2020
39dbb84
Merge branch 'master' into fix/shared/bounds-antimeridian
kodiakhq[bot] Jul 16, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
87 changes: 87 additions & 0 deletions packages/bathymetry/src/__test__/stac.test.ts
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,87 @@
import { GoogleTms } from '@basemaps/geo/build/tms/google';
import { LogConfig, LogType } from '@basemaps/shared';
import { round } from '@basemaps/test/build/rounding';
import o from 'ospec';
import { Hash } from '../hash';
import { Stac } from '../stac';

o.spec('stac', () => {
const origHash = Hash.hash;

o.afterEach(() => {
Hash.hash = origHash;
});
o('create', async () => {
(Hash as any).hash = (v: string): string => 'hash' + v;

const logger = LogConfig.get();

function name(a: string, b: string): string {
return a + '/' + b;
}

const bm = {
config: { tms: GoogleTms },
file: { name },
fileName: 'test-file.nc',
createSourceHash(l: LogType) {
return 'multihashResult' + (l === logger);
},
} as any;

const now = Date.now();

const stac = round((await Stac.create(bm, { x: 22, y: 33, z: 13 }, logger)) as any);

const date = Date.parse(stac.properties.datetime);

o(date >= now && date < now + 2000).equals(true);
o(/git@github\.com:linz\/basemaps.git#/.test(stac.links[1].href)).equals(true);
o(stac.links[1].rel).equals('derived_from');
o(/^\d+\.\d+\.\d+$/.test(stac.links[1].version)).equals(true);

o(stac).deepEquals({
// eslint-disable-next-line @typescript-eslint/camelcase
stac_version: '1.0.0',
// eslint-disable-next-line @typescript-eslint/camelcase
stac_extensions: ['proj', 'linz'],
id: '13-22-33',
type: 'Feature',
bounds: [-179.03320312, 84.92054529, -178.98925781, 84.92443459],
geometry: {
type: 'Polygon',
coordinates: [
[
[-179.03320312, 84.92054529],
[-179.03320312, 84.92443459],
[-178.98925781, 84.92443459],
[-178.98925781, 84.92054529],
[-179.03320312, 84.92054529],
],
],
},
properties: {
datetime: stac.properties.datetime,
'checksum:multihash': 'hashoutput/13-22-33',
'proj:epsg': 3857,
'linz:gdal:version': undefined,
},
assets: {
tiff: {
href: '13-22-33',
title: 'GeoTiff',
roles: ['data'],
type: 'image/vnd.stac.geotiff',
},
},
links: [
{
rel: 'derived_from',
href: 'test-file.nc',
'checksum:multihash': 'multihashResulttrue',
},
stac.links[1],
],
});
});
});
8 changes: 4 additions & 4 deletions packages/bathymetry/src/stac.ts
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
import { Tile, TileMatrixSet, GeoJson } from '@basemaps/geo';
import { Tile, TileMatrixSet } from '@basemaps/geo';
import { LogType, ProjectionTileMatrixSet } from '@basemaps/shared';
import * as cp from 'child_process';
import * as path from 'path';
Expand All @@ -20,8 +20,8 @@ async function create(bm: BathyMaker, tile: Tile, logger: LogType): Promise<Reco
const outputTiffPath = bm.file.name(FileType.Output, tileId);
const ptms = ProjectionTileMatrixSet.get(tms.projection.code);

const bounds = ptms.tileToWgs84Bounds(tile);
const geometry = GeoJson.toPolygon(bounds.toPolygon());
const bounds = ptms.tileToWgs84Bbox(tile);
const { geometry } = ptms.proj.boundsToGeoJsonFeature(ptms.tileToSourceBounds(tile));

const created = new Date().toISOString();
return {
Expand All @@ -31,7 +31,7 @@ async function create(bm: BathyMaker, tile: Tile, logger: LogType): Promise<Reco
stac_extensions: ['proj', 'linz'],
id: tileId,
type: 'Feature',
bounds: bounds.toBbox(),
bounds,
geometry,
properties: {
datetime: created,
Expand Down
70 changes: 70 additions & 0 deletions packages/shared/src/proj/__test__/projection.test.ts
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,76 @@ o.spec('Projection', () => {
]);
});

o.spec('boundsToGeoJsonFeature', () => {
o('simple', () => {
const ans = round(
googleProj.boundsToGeoJsonFeature(
{
x: -19929885.00696367,
y: 19871181.369240656,
width: 48921.969810251147,
height: 4891.969810251147,
},
{ name: '13-22-33' },
),
);

o(ans).deepEquals({
type: 'Feature',
geometry: {
type: 'Polygon',
coordinates: [
[
[-179.03320312, 84.92054529],
[-179.03320312, 84.92443459],
[-178.59372959, 84.92443459],
[-178.59372959, 84.92054529],
[-179.03320312, 84.92054529],
],
],
},
properties: { name: '13-22-33' },
});
});

o('crosses antimeridian', () => {
const ans = round(
nztmProj.boundsToGeoJsonFeature(
{ x: 1293760, y: 5412480, width: 1246880, height: 1146880 },
{ name: '1-2-1' },
),
);

o(ans).deepEquals({
type: 'Feature',
geometry: {
type: 'MultiPolygon',
coordinates: [
[
[
[169.33783769, -41.38093178],
[180, -41.22297022],
[180, -30.9080246],
[169.79094917, -31.05964823],
[169.33783769, -41.38093178],
],
],
[
[
[-175.8429206, -40.89543649],
[-177.19778597, -30.72635276],
[-180, -30.9080246],
[-180, -41.22297022],
[-175.8429206, -40.89543649],
],
],
],
},
properties: { name: '1-2-1' },
});
});
});

o('projectMultipolygon', () => {
const poly = [
Bounds.fromBbox([
Expand Down
56 changes: 7 additions & 49 deletions packages/shared/src/proj/__test__/projection.tile.matrix.set.test.ts
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -75,58 +75,16 @@ o.spec('ProjectionTileMatrixSet', () => {
});
});

o.spec('tileToWgs84Bounds', () => {
o('should convert base tiles', () => {
const pt = googlePtms.tileToWgs84Bounds({ x: 0, y: 0, z: 0 });
// GoogleTms extent is slightly off the epsg:3857 numbers
Approx.bounds(pt, Bounds.fromJson({ x: -180, y: -85.051, width: 360, height: 170.1 }), 'bounds', 0.01);

const lowerCorner = googlePtms.proj.toWsg84(GoogleTms.def.boundingBox.lowerCorner);
const upperCorner = googlePtms.proj.toWsg84(GoogleTms.def.boundingBox.upperCorner);
o.spec('tileToWgs84Bbox', () => {
o('should handle antimeridian', () => {
const pt = nztmPtms.tileToWgs84Bbox({ x: 2, y: 1, z: 1 });

Approx.equal(pt.x, lowerCorner[0], 'minX');
Approx.equal(pt.y, lowerCorner[1], 'minY');
Approx.equal(pt.right, upperCorner[0], 'maxX');
Approx.equal(pt.bottom, upperCorner[1], 'maxY');
o(round(pt)).deepEquals([170.05982382, -20.71836222, -179.34441046, -10.28396555]);
});

o('should round trip', () => {
const [minX, minY, maxX, maxY] = googlePtms.tileToWgs84Bounds({ x: 0, y: 0, z: 0 }).toBbox();
const topLeft = googlePtms.proj.fromWsg84([minX, maxY]);
const bottomRight = googlePtms.proj.fromWsg84([maxX, minY]);

// Lower as in number size, not lower as in position
Approx.equal(topLeft[0], GoogleTms.def.boundingBox.lowerCorner[0], 'topLeft.x');
Approx.equal(topLeft[1], GoogleTms.def.boundingBox.upperCorner[1], 'topLeft.y');

Approx.equal(bottomRight[0], GoogleTms.def.boundingBox.upperCorner[0], 'bottomRight.x');
Approx.equal(bottomRight[1], GoogleTms.def.boundingBox.lowerCorner[1], 'bottomRight.y');
});

o('should round trip bottom right individual tiles', () => {
const [minX, minY, maxX, maxY] = googlePtms.tileToWgs84Bounds({ x: 1, y: 1, z: 1 }).toBbox();
const topLeft = googlePtms.proj.fromWsg84([minX, maxY]);
const bottomRight = googlePtms.proj.fromWsg84([maxX, minY]);

// top left is in the center
Approx.equal(topLeft[0], 0, 'topLeft.x');
Approx.equal(topLeft[1], 0, 'topLeft.y');

Approx.equal(bottomRight[0], GoogleTms.def.boundingBox.upperCorner[0], 'bottomRight.x');
Approx.equal(bottomRight[1], GoogleTms.def.boundingBox.lowerCorner[1], 'bottomRight.y');
});

o('should round trip top left tile', () => {
const [minX, minY, maxX, maxY] = googlePtms.tileToWgs84Bounds({ x: 0, y: 0, z: 1 }).toBbox();
const topLeft = googlePtms.proj.fromWsg84([minX, maxY]);
const bottomRight = googlePtms.proj.fromWsg84([maxX, minY]);

// bottom right is in the center
Approx.equal(bottomRight[0], 0, 'bottomRight.x');
Approx.equal(bottomRight[1], 0, 'bottomRight.y');

Approx.equal(topLeft[0], GoogleTms.def.boundingBox.lowerCorner[0], 'topLeft.x');
Approx.equal(topLeft[1], GoogleTms.def.boundingBox.upperCorner[1], 'topLeft.y');
o('should convert base tiles', () => {
const pt = googlePtms.tileToWgs84Bbox({ x: 0, y: 0, z: 0 });
o(round(pt)).deepEquals([-180, -85.05112878, 180, 85.05112878]);
});
});

Expand Down
15 changes: 6 additions & 9 deletions packages/shared/src/proj/projection.tile.matrix.set.ts
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,17 +74,14 @@ export class ProjectionTileMatrixSet {
}

/** Convert a tile to the wgs84 bounds */
tileToWgs84Bounds(tile: Tile): Bounds {
tileToWgs84Bbox(tile: Tile): [number, number, number, number] {
const ul = this.tms.tileToSource(tile);
const lr = this.tms.tileToSource({ x: tile.x + 1, y: tile.y + 1, z: tile.z });
const [ulLon, ulLat] = this.proj.toWsg84([ul.x, ul.y]);
const [lrLon, lrLat] = this.proj.toWsg84([lr.x, lr.y]);
return new Bounds(
Math.min(ulLon, lrLon),
Math.min(ulLat, lrLat),
Math.abs(lrLon - ulLon),
Math.abs(lrLat - ulLat),
);

const [swLon, swLat] = this.proj.toWsg84([ul.x, lr.y]);
const [neLon, neLat] = this.proj.toWsg84([lr.x, ul.y]);

return [swLon, swLat, neLon, neLat];
}

/**
Expand Down
49 changes: 32 additions & 17 deletions packages/shared/src/proj/projection.ts
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
import { Bounds, Epsg, EpsgCode, GeoJson } from '@basemaps/geo';
import { BoundingBox, Epsg, EpsgCode, GeoJson } from '@basemaps/geo';
import { Position } from 'geojson';
import Proj from 'proj4';
import { NamedBounds } from '../aws/tile.metadata.base';
Expand Down Expand Up @@ -93,24 +93,39 @@ export class Projection {
return this.projection.inverse;
}

/** Convert a tile covering to a GeoJSON FeatureCollection */
toGeoJson(files: NamedBounds[]): GeoJSON.FeatureCollection {
/**
* Convert a source bounds to a WSG84 GeoJSON Feature

* @param bounds in source epsg
* @param properties any properties to include in the feature such as name
*/
boundsToGeoJsonFeature(bounds: BoundingBox, properties = {}): GeoJSON.Feature {
const { toWsg84 } = this;
const polygons: GeoJSON.Feature[] = [];
for (const bounds of files) {
const polygon = GeoJson.toFeaturePolygon(
[
Bounds.fromJson(bounds)
.toPolygon()[0]
.map((p) => toWsg84(p)),
],
{
name: bounds.name,
},
);
polygons.push(polygon);
const sw = toWsg84([bounds.x, bounds.y]);
const se = toWsg84([bounds.x + bounds.width, bounds.y]);
const nw = toWsg84([bounds.x, bounds.y + bounds.height]);
const ne = toWsg84([bounds.x + bounds.width, bounds.y + bounds.height]);

if (sw[0] < se[0] && nw[0] < ne[0]) {
return GeoJson.toFeaturePolygon([[sw, nw, ne, se, sw]], properties);
}

return GeoJson.toFeatureCollection(polygons);
// crosses antimeridian so need to split polygon at antimeridian

// calculate where antimeridian is at north and south bounds
const northFraction = (180 - nw[0]) / (ne[0] + 360 - nw[0]);
const southFraction = (180 - sw[0]) / (se[0] + 360 - sw[0]);
const n180 = toWsg84([bounds.x + bounds.width * northFraction, bounds.y + bounds.height])[1];
const s180 = toWsg84([bounds.x + bounds.width * southFraction, bounds.y])[1];

return GeoJson.toFeatureMultiPolygon(
[[[sw, [180, s180], [180, n180], nw, sw]], [[se, ne, [-180, n180], [-180, s180], se]]],
properties,
);
}

/** Convert a tile covering to a GeoJSON FeatureCollection */
toGeoJson(files: NamedBounds[]): GeoJSON.FeatureCollection {
return GeoJson.toFeatureCollection(files.map((f) => this.boundsToGeoJsonFeature(f, { name: f.name })));
}
}