feat(geojson): Improve GeoJSON compliance

packages/bathymetry/src/stac.ts

@@ -21,7 +21,7 @@ async function create(bm: BathyMaker, tile: Tile, logger: LogType): Promise<Reco
     const ptms = ProjectionTileMatrixSet.get(tms.projection.code);
 
     const bounds = ptms.tileToWgs84Bbox(tile);
-    const { geometry } = ptms.proj.boundsToGeoJsonFeature(ptms.tileToSourceBounds(tile));
+    const { geometry } = ptms.proj.boundsToGeoJsonFeature(tms.tileToSourceBounds(tile));
 
     const created = new Date().toISOString();
     return {

packages/cli/package.json

@@ -26,6 +26,7 @@
     }
   ],
   "dependencies": {
+    "@linzjs/geojson": "^1.0.0",
     "@basemaps/geo": "^4.4.0",
     "@basemaps/shared": "^4.5.0",
     "@cogeotiff/core": "^2.2.0",
@@ -34,9 +35,7 @@
     "@cogeotiff/source-url": "^2.2.0",
     "@rushstack/ts-command-line": "^4.3.13",
     "ansi-colors": "^4.1.1",
-    "lineclip": "^1.1.5",
     "p-limit": "^3.0.1",
-    "polygon-clipping": "^0.15.1",
     "pretty-json-log": "^0.3.1"
   },
   "devDependencies": {

packages/cli/src/cog/__test__/cog.test.ts

@@ -33,7 +33,7 @@ o.spec('cog', () => {
 
             const name = '4-15-10';
 
-            job.files = [{ name, ...targetPtms.tileToSourceBounds({ x: 15, y: 10, z: 4 }) }];
+            job.files = [{ name, ...targetPtms.tms.tileToSourceBounds({ x: 15, y: 10, z: 4 }) }];
 
             await buildCogForName(job, name, '/tmp/test.vrt', '/tmp/out-tiff', logger, true);
             o(convertArgs[0].info).equals(logger.info);

packages/cli/src/cog/__test__/cutline.test.ts

@@ -1,11 +1,12 @@
-import { EpsgCode } from '@basemaps/geo';
+import { EpsgCode, Bounds } from '@basemaps/geo';
 import { ProjectionTileMatrixSet } from '@basemaps/shared';
 import { qkToName } from '@basemaps/shared/build/proj/__test__/test.util';
 import { round } from '@basemaps/test/build/rounding';
 import o from 'ospec';
-import { Cutline } from '../cutline';
+import { Cutline, polyContainsBounds } from '../cutline';
 import { SourceMetadata } from '../types';
 import { SourceTiffTestHelper } from './source.tiff.testhelper';
+import { MultiPolygon } from '@linzjs/geojson';
 
 o.spec('cutline', () => {
     const testDir = `${__dirname}/../../../__test.assets__`;
@@ -32,6 +33,47 @@ o.spec('cutline', () => {
         });
     });
 
+    o('polyContainsBounds', () => {
+        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(polys, Bounds.fromBbox([-6, -6, -3, -3]))).equals(true);
+
+        o(polyContainsBounds(polys, Bounds.fromBbox([-3, -4, 4, 4]))).equals(false);
+        o(polyContainsBounds(polys, Bounds.fromBbox([6, -8, 5, 5]))).equals(false);
+    });
+
     o('loadCutline', async () => {
         const cutline = new Cutline(googlePtms, await Cutline.loadCutline(testDir + '/mana.geojson'));
         const geojson = round(cutline.toGeoJson());

packages/cli/src/cog/cutline.ts

@@ -1,12 +1,10 @@
 import { Bounds, Epsg, GeoJson, Tile, TileMatrixSet } from '@basemaps/geo';
 import { compareName, FileOperator, NamedBounds, ProjectionTileMatrixSet } from '@basemaps/shared';
+import { Projection } from '@basemaps/shared/build/proj/projection';
+import { clipMultipolygon, intersection, MultiPolygon, Ring, union } from '@linzjs/geojson';
 import { FeatureCollection } from 'geojson';
 import { CoveringFraction, MaxImagePixelWidth } from './constants';
 import { CogJob, SourceMetadata } from './types';
-import { clipMultipolygon, polyContainsBounds } from './clipped.multipolygon';
-import pc, { MultiPolygon, Ring, Polygon } from 'polygon-clipping';
-import { Projection } from '@basemaps/shared/build/proj/projection';
-const { intersection, union } = pc;
 
 export interface FeatureCollectionWithCrs extends FeatureCollection {
     crs: {
@@ -33,6 +31,13 @@ function namedBounds(tms: TileMatrixSet, tile: Tile): NamedBounds {
     return { name: TileMatrixSet.tileToName(tile), ...tms.tileToSourceBounds(tile).toJson() };
 }
 
+export function polyContainsBounds(poly: MultiPolygon, bounds: Bounds): boolean {
+    const clipped = clipMultipolygon(poly, bounds.toBbox());
+    if (clipped.length != 1 || clipped[0].length != 1 || clipped[0][0].length != 5) return false;
+
+    return Bounds.fromMultiPolygon(clipped).containsBounds(bounds);
+}
+
 /**
  * Filter out duplicate tiles
  */
@@ -53,12 +58,16 @@ function addNonDupes(list: Tile[], addList: Tile[]): void {
 }
 
 export class Cutline {
+    /** The polygon to clip source imagery to */
     clipPoly: MultiPolygon = [];
     targetPtms: ProjectionTileMatrixSet;
+    tms: TileMatrixSet; // convience to targetPtms.tms
+    /** How much blending to apply at the clip line boundary */
     blend: number;
+    /** For just one cog to cover the imagery */
     oneCog: boolean;
-    tms: TileMatrixSet; // convience to targetPtms.tms
-    private srcPoly: MultiPolygon = [];
+    /** the polygon outlining a area covered by the source imagery and clip polygon */
+    srcPoly: MultiPolygon = [];
 
     /**
      * Create a Cutline instance from a `GeoJSON FeatureCollection`.
@@ -128,7 +137,7 @@ export class Cutline {
         const tile = TileMatrixSet.nameToTile(name);
         const sourceCode = Projection.get(job.source.projection);
         const targetCode = this.targetPtms.proj;
-        const tileBounds = this.targetPtms.tileToSourceBounds(tile);
+        const tileBounds = this.tms.tileToSourceBounds(tile);
         const tilePadded = this.padBounds(tileBounds, job.source.resZoom);
 
         let tileBoundsInSrcProj = tilePadded;
@@ -139,8 +148,9 @@ export class Cutline {
             tileBoundsInSrcProj = Bounds.fromMultiPolygon(poly);
         }
 
+        const paddedBbox = tilePadded.toBbox();
         if (this.clipPoly.length > 0) {
-            const poly = clipMultipolygon(this.clipPoly, tilePadded);
+            const poly = clipMultipolygon(this.clipPoly, paddedBbox);
             if (poly.length == 0) {
                 // this tile is not needed
                 this.clipPoly = [];
@@ -231,7 +241,7 @@ export class Cutline {
         minZ: number,
         coveringFraction: number,
     ): { tiles: Tile[]; fractionCovered: number } {
-        const clipBounds = this.targetPtms.tileToSourceBounds(tile);
+        const clipBounds = this.tms.tileToSourceBounds(tile).toBbox();
 
         srcArea = clipMultipolygon(srcArea, clipBounds);
 
@@ -279,7 +289,7 @@ export class Cutline {
 
         // merge imagery bounds
         for (const image of sourceMetadata.bounds) {
-            const poly = Bounds.fromJson(image).scaleFromCenter(SourceSmoothScale).toPolygon() as Polygon;
+            const poly = [Bounds.fromJson(image).scaleFromCenter(SourceSmoothScale).toPolygon()] as MultiPolygon;
             srcPoly = union(srcPoly, poly);
         }
 
@@ -294,7 +304,7 @@ export class Cutline {
         if (this.clipPoly.length == 0) return;
 
         const srcBounds = Bounds.fromMultiPolygon(srcPoly);
-        const boundsPadded = this.padBounds(srcBounds, resZoom);
+        const boundsPadded = this.padBounds(srcBounds, resZoom).toBbox();
 
         const poly = clipMultipolygon(this.clipPoly, boundsPadded);
         if (poly.length == 0) {

packages/cli/src/gdal/gdal.config.ts

@@ -1,4 +1,5 @@
 import { Epsg } from '@basemaps/geo';
+import { BBox } from '@linzjs/geojson';
 
 export type GdalCogBuilderOptionsResampling =
     | 'nearest'
@@ -41,7 +42,7 @@ export interface GdalCogBuilderOptions {
 
     /** Limit the output to a bounding box
      */
-    bbox?: [number, number, number, number];
+    bbox?: BBox;
 
     /**
      * Compression to use for the cog

packages/geo/src/bounds.ts

@@ -174,7 +174,8 @@ export class Bounds implements BoundingBox {
     }
 
     /**
-     * Find the bounds of a GeoJson MultiPolygon
+     * Find the bounds of a MultiPolygon.
+     * Does not work with WGS84 when crosses antimeridian.
 
      * @param multipoly the polygon to measure
      */
@@ -198,7 +199,7 @@ export class Bounds implements BoundingBox {
 
     /**
      * Convert a BBox(minX, minY, maxX, maxY) to Bounds(x,y, width, height).
-     * Takes into account the antimeridian.
+     * Does not work with WGS84 when crosses the antimeridian.
      * @param bbox
      */
     public static fromBbox([x1, y1, x2, y2]: number[]): Bounds {

packages/geo/src/index.ts

@@ -1,6 +1,6 @@
-export { Bounds, BoundingBox, Size } from './bounds';
-export { QuadKey } from './quad.key';
-export { GeoJson } from './geo.json';
+export { BoundingBox, Bounds, Size } from './bounds';
 export { Epsg, EpsgCode } from './epsg';
-export { TileMatrixSet, Tile } from './tile.matrix.set';
+export { GeoJson } from './geo.json';
+export { QuadKey } from './quad.key';
+export { Tile, TileMatrixSet } from './tile.matrix.set';
 export { WmtsLayer, WmtsProvider } from './wmts/wmts';

packages/lambda-xyz/package.json

@@ -6,6 +6,7 @@
   "author": "",
   "license": "MIT",
   "dependencies": {
+    "@linzjs/geojson": "^1.0.0",
     "@basemaps/geo": "^4.4.0",
     "@basemaps/lambda": "^4.5.0",
     "@basemaps/shared": "^4.5.0",

packages/lambda-xyz/src/__test__/wmts.capability.test.ts

@@ -108,7 +108,7 @@ o.spec('WmtsCapabilities', () => {
         o(xml).deepEquals('<?xml version="1.0"?>\n' + raw?.toString());
 
         o(createHash('sha256').update(Buffer.from(xml)).digest('base64')).equals(
-            'LeUIcNjN/lFznAoFvrDzOh9uBbxC7nn+bcbKD9FogoA=',
+            'qBC5TNmXDmxYhIS2SAJQwwWrHxXKjm7sk5K29LAOtt8=',
         );
     });
 

packages/lambda-xyz/src/__test__/xyz.test.ts

@@ -157,11 +157,15 @@ o.spec('LambdaXyz', () => {
         });
 
         o('should 304 if a xml is not modified', async () => {
-            const key = 'biarWiiP+sp+4QsJnuQwxlxW3zEnipGptLywav1E7Cs=';
+            const key = 'oxWjinmkGeDsEoFBW1wZ1cUTXD1yth4gkJp5EsphoU8=';
             const request = mockRequest('/v1/tiles/WMTSCapabilities.xml', 'get', { 'if-none-match': key });
 
             const res = await handleRequest(request);
-            if (res.status == 200) o(res.header('eTaG')).equals(key); // this line is useful for discovering the new etag
+            if (res.status == 200) {
+                o(res.header('eTaG')).equals(key); // this line is useful for discovering the new etag
+                return;
+            }
+
             o(res.status).equals(304);
             o(rasterMock.calls.length).equals(0);
 

packages/lambda-xyz/src/wmts.capability.ts

@@ -1,10 +1,11 @@
-import { Epsg, TileMatrixSet, WmtsLayer, WmtsProvider, Bounds } from '@basemaps/geo';
+import { Bounds, Epsg, TileMatrixSet, WmtsLayer, WmtsProvider } from '@basemaps/geo';
 import { GoogleTms } from '@basemaps/geo/build/tms/google';
 import { Nztm2000Tms } from '@basemaps/geo/build/tms/nztm2000';
 import { TileMetadataProviderRecord, V, VNodeElement } from '@basemaps/shared';
+import { Projection } from '@basemaps/shared/build/proj/projection';
 import { ImageFormatOrder } from '@basemaps/tiler';
+import { BBox, Wgs84 } from '@linzjs/geojson';
 import { TileSet } from './tile.set';
-import { Projection } from '@basemaps/shared/build/proj/projection';
 
 function getTileMatrixSet(projection: Epsg): TileMatrixSet {
     switch (projection) {
@@ -41,8 +42,8 @@ const CapabilitiesAttrs = {
     version: '1.0.0',
 };
 
-function wgs84Extent(layer: WmtsLayer): Bounds {
-    return Projection.get(layer.projection.code).boundsToWgs84(layer.extent);
+function wgs84Extent(layer: WmtsLayer): BBox {
+    return Projection.get(layer.projection.code).boundsToWgs84BoundingBox(layer.extent);
 }
 
 export class WmtsCapabilities {
@@ -72,12 +73,11 @@ export class WmtsCapabilities {
     }
 
     buildWgs84BoundingBox(layers: WmtsLayer[], tagName = 'ows:WGS84BoundingBox'): VNodeElement {
-        let bounds = wgs84Extent(layers[0]);
+        let bbox = wgs84Extent(layers[0]);
         for (let i = 1; i < layers.length; ++i) {
-            bounds = bounds.union(wgs84Extent(layers[i]));
+            bbox = Wgs84.union(bbox, wgs84Extent(layers[i]));
         }
 
-        const bbox = bounds.toBbox();
         return V(
             tagName,
             { crs: 'urn:ogc:def:crs:OGC:2:84' },

packages/linzjs-geojson/CHANGELOG.md

@@ -0,0 +1,4 @@
+# Change Log
+
+All notable changes to this project will be documented in this file.
+See [Conventional Commits](https://conventionalcommits.org) for commit guidelines.

packages/linzjs-geojson/README.md

@@ -0,0 +1,34 @@
+# GeoJSON computation
+
+Utility functions for working with GeoJSON multi polygons and bounding boxes. In particular for handling the anti-meridian.
+
+## Usage
+
+
+### Wgs84
+```javascript
+import { Wgs84 } from '@linzjs/geojson';
+
+assert(Wgs84.normLon(-163.12345 - 720) == -163.12345;
+
+assert(Wgs84.crossesAM(-175, 175));
+
+assert(Wgs84.delta(-175, 170) == -15);
+
+assert(deepEqual(Wgs84.union([175, -42, -178, -41], [-170, -43, -160, -42]), [175, -43, -160, -41]));
+```
+
+### MultiPolygon
+
+```javascript
+import { clipMultipolygon,  multiPolygonToWgs84 } from '@linzjs/geojson';
+import Proj from 'proj4';
+
+// polygons clipped to bounding box; no degenerate edges
+const clipped = clipMultipolygon(polygons, [-2, -2, 1, 1]);
+
+const nztmToWgs84 = Proj('epsg:2193', 'epsg:4326').forward;
+
+// nztm polygons converted to wgs84 split at anti-meridian
+const splitPolys = multiPolygonToWgs84(nztmPolygons, nztmToWgs84);
+```