From 6c832607e137f2a740c1f1960b9e7c7c66a5c2bf Mon Sep 17 00:00:00 2001 From: Geoff Jacobsen Date: Wed, 15 Jul 2020 11:58:51 +1200 Subject: [PATCH] fix(shared): handle bounds crossing antimeridian Ensure compilance with GeoJSON RFC7946: 1. Correct order of Bounding box coordinates 2. Split bounding box poygon in two when crossing antimeridian --- packages/bathymetry/src/__test__/stac.test.ts | 87 +++++++++++++++++++ packages/bathymetry/src/stac.ts | 8 +- .../src/proj/__test__/projection.test.ts | 70 +++++++++++++++ .../projection.tile.matrix.set.test.ts | 56 ++---------- .../src/proj/projection.tile.matrix.set.ts | 15 ++-- packages/shared/src/proj/projection.ts | 49 +++++++---- 6 files changed, 206 insertions(+), 79 deletions(-) create mode 100644 packages/bathymetry/src/__test__/stac.test.ts diff --git a/packages/bathymetry/src/__test__/stac.test.ts b/packages/bathymetry/src/__test__/stac.test.ts new file mode 100644 index 000000000..613cb58c7 --- /dev/null +++ b/packages/bathymetry/src/__test__/stac.test.ts @@ -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], + ], + }); + }); +}); diff --git a/packages/bathymetry/src/stac.ts b/packages/bathymetry/src/stac.ts index c08137dd2..e2663ca89 100644 --- a/packages/bathymetry/src/stac.ts +++ b/packages/bathymetry/src/stac.ts @@ -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'; @@ -20,8 +20,8 @@ async function create(bm: BathyMaker, tile: Tile, logger: LogType): Promise { ]); }); + 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([ diff --git a/packages/shared/src/proj/__test__/projection.tile.matrix.set.test.ts b/packages/shared/src/proj/__test__/projection.tile.matrix.set.test.ts index 1cdc6af25..8ea5a6a68 100644 --- a/packages/shared/src/proj/__test__/projection.tile.matrix.set.test.ts +++ b/packages/shared/src/proj/__test__/projection.tile.matrix.set.test.ts @@ -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]); }); }); diff --git a/packages/shared/src/proj/projection.tile.matrix.set.ts b/packages/shared/src/proj/projection.tile.matrix.set.ts index 5f728518c..caaca8002 100644 --- a/packages/shared/src/proj/projection.tile.matrix.set.ts +++ b/packages/shared/src/proj/projection.tile.matrix.set.ts @@ -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]; } /** diff --git a/packages/shared/src/proj/projection.ts b/packages/shared/src/proj/projection.ts index 8d8f9ae42..2b2047c0f 100644 --- a/packages/shared/src/proj/projection.ts +++ b/packages/shared/src/proj/projection.ts @@ -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'; @@ -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 }))); } }