feat: attempt to remove slivers by dropping small polygons TDE-1130 (#…

…1119)

#### Motivation

Quite a few of our datasets are aligned to a EPSG:2193 tile grid, when
reprojecting these area polygons to a EPSG:4326, small slithers can
appear

for example in this dataset where one dataset (magenta) has 5 tiles
across and one has 4 (blue) only the two corner points are stored

![image
(41)](https://github.com/user-attachments/assets/c3302f71-52af-473e-9c77-e08bf481d0ed)

As you zoom in on the corner slithers appear

![image
(42)](https://github.com/user-attachments/assets/0202fecd-b564-4056-9421-17e55d9a4e3f)

#### Modification

Find polygons that are very small and drop them.

#### Checklist

_If not applicable, provide explanation of why._

- [ ] Tests updated
- [ ] Docs updated
- [ ] Issue linked in Title

package.json

@@ -46,6 +46,7 @@
     "@linzjs/geojson": "^7.10.0",
     "@octokit/core": "^5.0.0",
     "@octokit/plugin-rest-endpoint-methods": "^10.1.1",
+    "@turf/buffer": "^7.1.0",
     "@wtrpc/core": "^1.0.2",
     "ajv": "^8.17.1",
     "ajv-formats": "^2.1.1",

src/commands/mapsheet-coverage/__test__/mapsheet.coverage.test.ts

@@ -8,7 +8,7 @@ import { featuresToMultiPolygon } from '@linzjs/geojson';
 import { StacCollection } from 'stac-ts';
 
 import { MapSheet } from '../../../utils/mapsheet.js';
-import { commandMapSheetCoverage } from '../mapsheet.coverage.js';
+import { commandMapSheetCoverage, isLargeRegion } from '../mapsheet.coverage.js';
 
 // convert a collection of map sheets into a multipolygon
 function mapSheetToGeoJson(...sheetCodes: string[]): GeoJSON.Feature {
@@ -200,3 +200,72 @@ describe('mapsheet-coverage', () => {
     assert.equal(captureDates.features.length, 1);
   });
 });
+
+describe('isLargeRegion', () => {
+  const OneMetreInDegrees = 1e-5;
+
+  it('should not remove large polygons', () => {
+    // this polygon is 1 degree x 1 degree at approx 110KM x 110KM
+    assert.equal(
+      isLargeRegion([
+        [
+          [0, 0],
+          [0, -1],
+          [-1, -1],
+          [-1, 0],
+          [0, 0],
+        ],
+      ]),
+      true,
+    );
+  });
+
+  it('should remove polygons with long slivers', () => {
+    // this polygon is 1 degree x 1cm very small slivers
+    // which will be destroyed by the buffering inwards
+    assert.equal(
+      isLargeRegion([
+        [
+          [0, 0],
+          [1, 0],
+          [1, OneMetreInDegrees * 0.01],
+          [0, OneMetreInDegrees * 0.01],
+          [0, 0],
+        ],
+      ]),
+      false,
+    );
+  });
+
+  it('should remove small squares', () => {
+    // this polygon is 1m x 1m, with a area of 1E-10, it should be removed
+    assert.equal(
+      isLargeRegion([
+        [
+          [0, 0],
+          [OneMetreInDegrees, 0],
+          [OneMetreInDegrees, OneMetreInDegrees],
+          [0, OneMetreInDegrees],
+          [0, 0],
+        ],
+      ]),
+      false,
+    );
+  });
+
+  it('should keep medium sized polygons', () => {
+    // this polygon is 100m x 100m
+    assert.equal(
+      isLargeRegion([
+        [
+          [0, 0],
+          [100 * OneMetreInDegrees, 0],
+          [100 * OneMetreInDegrees, 100 * OneMetreInDegrees],
+          [0, 100 * OneMetreInDegrees],
+          [0, 0],
+        ],
+      ]),
+      true,
+    );
+  });
+});

src/commands/mapsheet-coverage/mapsheet.coverage.ts

@@ -3,7 +3,8 @@ import { gzipSync } from 'node:zlib';
 import { ConfigTileSetRaster } from '@basemaps/config';
 import { EpsgCode } from '@basemaps/geo';
 import { fsa } from '@chunkd/fs';
-import { truncate } from '@linzjs/geojson';
+import { Area, truncate } from '@linzjs/geojson';
+import buffer from '@turf/buffer';
 import { command, number, option, optional, string } from 'cmd-ts';
 import pLimit from 'p-limit';
 import { basename } from 'path/posix';
@@ -41,6 +42,39 @@ function forceMultiPolygon(f: GeoJSON.Feature): GeoJSON.Feature<GeoJSON.MultiPol
   throw new Error(`${f.geometry.type}: is not a polygon`);
 }
 
+/**
+ * a polygon of 0.00001 x 0.00001 (about 1m x 1m) converts to 1e-10 degrees^2
+ *
+ * @warning
+ * This is a approximation and should not be used for anything other than rough scale comparisons
+ */
+const DegreesSquareToMetersSquared = 1e-10;
+
+/**
+ * Determine if the polygon is a large(ish) region,
+ *
+ * @param poly Polygon in EPSG:4326 (WGS84 lat, lon)
+ * @param areaInMetersSquared Polygons that are at least this big are always included without further testing
+ *
+ * @returns if the polygon is considered large
+ */
+export function isLargeRegion(poly: pc.Polygon, areaInMetersSquared: number = 1_000): boolean {
+  const polyArea = Area.polygon(poly);
+  const smallArea = areaInMetersSquared * DegreesSquareToMetersSquared;
+
+  // Area is stored as square degrees, 1e-10 is approx 1m^2 (?)
+  // TODO: is there a better number than 1e-6, could scale it from the GSD of the dataset
+  if (polyArea > smallArea) return true;
+
+  // If the buffer in destroys the polygon this polygon likely isn't really useful
+  // TODO: 1m works when the dataset is 1m, but not all datasets are 1m resolution
+  const bufferedIn = buffer({ type: 'Polygon', coordinates: poly }, -1, { units: 'meters' });
+  if (bufferedIn == null) return false;
+
+  // TODO: should the buffered dataset be used
+  return true;
+}
+
 export const commandMapSheetCoverage = command({
   name: 'mapsheet-coverage',
   description: 'Create a list of mapsheets needing to be created from a basemaps configuration',
@@ -106,7 +140,7 @@ export const commandMapSheetCoverage = command({
     const captureDates = { type: 'FeatureCollection', features: [] as GeoJSON.Feature[] };
 
     // All the coordinates currently used
-    let layersCombined: pc.MultiPolygon[] = [];
+    let layersCombined: pc.MultiPolygon = [];
 
     // MapSheetName to List of source files required
     const mapSheets = new Map<string, string[]>();
@@ -175,10 +209,13 @@ export const commandMapSheetCoverage = command({
       );
 
       // GeoJSON over about 8 decimal places starts running into floating point math errors
-      truncate(captureArea);
+      truncate(captureArea, 8);
 
+      const diff = [];
       // Determine if this layer has any additional information to the existing layers
-      const diff = pc.difference(captureArea.geometry.coordinates as pc.MultiPolygon, layersCombined);
+      for (const poly of pc.difference(captureArea.geometry.coordinates as pc.MultiPolygon, layersCombined)) {
+        if (isLargeRegion(poly)) diff.push(poly);
+      }
 
       // Layer has no information included in the output
       if (diff.length === 0) {
@@ -208,6 +245,7 @@ export const commandMapSheetCoverage = command({
       }
 
       layersCombined = pc.union(layersCombined, captureArea.geometry.coordinates as pc.MultiPolygon);
+
       captureDates.features.push({
         type: 'Feature',
         geometry: { type: 'MultiPolygon', coordinates: diff },