Skip to content

Commit

Permalink
fix(shared): handle bounds crossing antimeridian
Browse files Browse the repository at this point in the history 
Ensure compilance with GeoJSON RFC7946:

 1. Correct order of Bounding box coordinates
 2. Split bounding box poygon in two when crossing antimeridian
  • Loading branch information
Geoff Jacobsen committed Jul 16, 2020
1 parent 11819bb commit 6c83260
Show file tree
Hide file tree
Showing 6 changed files with 206 additions and 79 deletions.
87 changes: 87 additions & 0 deletions packages/bathymetry/src/__test__/stac.test.ts
View file
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
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
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
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
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
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 })));
}
}

0 comments on commit 6c83260

Please sign in to comment.