Preprocessor update (#393)

* add missing toJsonFile helper

* getBiggestPolygon: break out into helper function

* change getBiggestPolygon to biggestPolygon

* polygonSmokeTest simplify thrown error message

* polygonSmokeTest: only log on ValidationError

* polygonSmokeTest: tune logging

* update preprocessing functions to use new global datasources (needs to be changed at end)

* switch back to production global-datasources urls, update e2e tests

* clipToGeography - increase test timeout

packages/base-project/src/util/clipToGeography.test.ts

@@ -480,7 +480,7 @@ describe("clipToGeography", () => {
     const clippedSketchBox = clippedSketch.bbox || bbox(clippedSketch);
     expect(clippedSketchArea).toEqual(sketchArea);
     expect(sketchBox).toEqual(clippedSketchBox);
-  });
+  }, 10_000);
 
   test("clipToGeography - no overlap", async () => {
     const curGeography = project.getGeographyById("world");

packages/geoprocessing/scripts/base/datasources/precalcRasterDatasource.test.e2e.ts

@@ -406,7 +406,7 @@ describe("precalcRasterDatasource", () => {
       metrics,
       (m) => m.geographyId === "geog-box-filter" && m.metricId === "sum",
     );
-    expect(boxFilterMetric.value).toEqual(70);
+    expect(boxFilterMetric.value).toEqual(68);
 
     const singleFilterMetric = firstMatchingMetric(
       metrics,
@@ -418,7 +418,7 @@ describe("precalcRasterDatasource", () => {
       metrics,
       (m) => m.geographyId === "geog-double-filter" && m.metricId === "sum",
     );
-    expect(doubleFilterMetric.value).toEqual(69);
+    expect(doubleFilterMetric.value).toEqual(67);
 
     fs.removeSync(dsFilePath);
     fs.removeSync(path.join(dstPath, `${rasterDatasourceId}.tif`));

packages/geoprocessing/scripts/base/datasources/precalcVectorDatasource.test.e2e.ts

@@ -818,21 +818,21 @@ describe("precalcDatasources", () => {
         (m) => m.geographyId === "geog-box-filter" && m.metricId === "area",
       );
       // Largest area value
-      expect(boxFilterMetric.value).toEqual(59_556_498_695.328_66);
+      expect(boxFilterMetric.value).toEqual(51_247_161_473.9979);
 
       const singleFilterMetric = firstMatchingMetric(
         metrics,
         (m) => m.geographyId === "geog-single-filter" && m.metricId === "area",
       );
       // Smallest area value, samoa only
-      expect(singleFilterMetric.value).toEqual(35_442_309_711.542_16);
+      expect(singleFilterMetric.value).toEqual(35_442_309_366.529_57);
 
       const doubleFilterMetric = firstMatchingMetric(
         metrics,
         (m) => m.geographyId === "geog-double-filter" && m.metricId === "area",
       );
       // Slightly larger area value, both samoa
-      expect(doubleFilterMetric.value).toEqual(37_350_139_043.666_25);
+      expect(doubleFilterMetric.value).toEqual(37_350_103_043.708_23);
 
       fs.removeSync(dsFilePath);
       fs.removeSync(path.join(dstPath, `${internalDatasourceId}.fgb`));

packages/geoprocessing/scripts/testing/polygonSmokeTest.ts

@@ -55,9 +55,14 @@ export const polygonSmokeTest = (
               example?.properties?.name,
             );
           } catch (error) {
-            console.log("error", example?.properties?.name, error);
+            console.log(
+              "preprocessor",
+              preprocessorName,
+              example?.properties?.name,
+            );
             if (error instanceof ValidationError) {
-              // ValidationErrors don't indicate failures, just comprehensive tests
+              console.log(error.message);
+              // ValidationErrors don't indicate failures, just comprehensive tests, so swallow and don't generate output
             } else {
               throw error;
             }

packages/geoprocessing/src/dataproviders/flatgeobuf.ts

@@ -28,7 +28,7 @@ export function fgBoundingBox(box: BBox): FgBoundingBox {
  * @param url url of flatgeobuf file
  * @param bbox optional bounding box to fetch features that intersect with
  * @returns feature array
- * @deprecated Use `loadCog` instead.
+ * @deprecated Use `loadFgb` instead.
  */
 export async function fgbFetchAll<F extends Feature<Geometry>>(
   url: string,

packages/geoprocessing/src/helpers/biggestPolygon.ts

@@ -0,0 +1,22 @@
+import { area, flatten } from "@turf/turf";
+import { Feature, MultiPolygon, Polygon } from "geojson";
+
+/**
+ * If feature is a MultiPolygon, scans and returns the polygon with the largest area.
+ */
+export function biggestPolygon(feature: Feature<Polygon | MultiPolygon>) {
+  if (feature.geometry.type === "MultiPolygon") {
+    // If multipolygon clip result, keep only the biggest piece
+    const flattened = flatten(feature);
+    let biggest = [0, 0];
+    for (let i = 0; i < flattened.features.length; i++) {
+      const a = area(flattened.features[i]);
+      if (a > biggest[0]) {
+        biggest = [a, i];
+      }
+    }
+    return flattened.features[biggest[1]] as Feature<Polygon>;
+  } else {
+    return feature;
+  }
+}

packages/geoprocessing/src/helpers/index.ts

@@ -13,6 +13,7 @@ export * from "./randomData.js";
 export * from "./sketch.js";
 export * from "./string.js";
 export * from "./ts.js";
+export * from "./fs.js";
 export * from "./units.js";
 export * from "./valueFormatter.js";
 export * from "./extraParams.js";
@@ -24,3 +25,4 @@ export * from "./callWithRetry.js";
 export * from "./toMultiPolygon.js";
 export * from "./removeOverlap.js";
 export * from "./removeHoles.js";
+export * from "./biggestPolygon.js";

packages/template-blank-project/src/functions/clipToLand.test.ts

@@ -1,6 +1,6 @@
 import { describe, test, expect, vi, afterEach, assert } from "vitest";
 import { clipToLand } from "./clipToLand.js";
-import { area, featureCollection, polygon } from "@turf/turf";
+import { area, polygon } from "@turf/turf";
 
 const landFeature = polygon([
   [
@@ -50,23 +50,20 @@ describe("clipToLand", () => {
     vi.restoreAllMocks();
   });
 
-  // Mock VectorDataSource fetchUnion method to return landFeature
+  // Mock loadFgb method to return landFeature
   // @ts-ignore
   vi.mock(import("@seasketch/geoprocessing"), async (importOriginal) => {
     const actual = await importOriginal();
-    const VectorDataSource = vi.fn();
-    VectorDataSource.prototype.fetchUnion = vi.fn(() =>
-      featureCollection([landFeature]),
-    );
-    return { ...actual, VectorDataSource };
+    const loadFgb = vi.fn(() => [landFeature]);
+    return { ...actual, loadFgb };
   });
 
   test("clipToLand - feature outside of land should throw", async () => {
     try {
       await clipToLand(outside);
     } catch (error: unknown) {
       if (error instanceof Error) {
-        expect(error.message).toBe("Feature is outside of boundary");
+        expect(error.message).toBe("Feature is outside of land boundary");
         return;
       }
     }

packages/template-blank-project/src/functions/clipToLand.ts

@@ -4,11 +4,14 @@ import {
   Sketch,
   isPolygonFeature,
   ValidationError,
-  VectorDataSource,
-  FeatureClipOperation,
-  clipToPolygonFeatures,
+  Polygon,
+  MultiPolygon,
+  loadFgb,
+  clipMultiMerge,
+  ensureValidPolygon,
+  biggestPolygon,
 } from "@seasketch/geoprocessing";
-import { bbox } from "@turf/turf";
+import { area, bbox, featureCollection } from "@turf/turf";
 
 /**
  * Preprocessor takes a Polygon feature/sketch and returns the portion that
@@ -18,24 +21,46 @@ export async function clipToLand(feature: Feature | Sketch): Promise<Feature> {
   if (!isPolygonFeature(feature)) {
     throw new ValidationError("Input must be a polygon");
   }
+
+  // throws if not valid with specific message
+  ensureValidPolygon(feature, {
+    minSize: 1,
+    enforceMinSize: false,
+    maxSize: 500_000 * 1000 ** 2, // Default 500,000 KM
+    enforceMaxSize: false,
+  });
+
   const featureBox = bbox(feature);
 
-  // Get land polygons - osm land vector datasource
-  const landDatasource = new VectorDataSource(
-    "https://d3p1dsef9f0gjr.cloudfront.net/",
+  // Get features from land datasource
+
+  const landFeatures: Feature<Polygon | MultiPolygon>[] = await loadFgb(
+    "https://gp-global-datasources-datasets.s3.us-west-1.amazonaws.com/global-coastline-daylight-v158.fgb",
+    featureBox,
   );
-  // one gid assigned per country, use to union subdivided pieces back together on fetch, prevents slivers
-  const landFC = await landDatasource.fetchUnion(featureBox, "gid");
 
-  const keepLand: FeatureClipOperation = {
-    operation: "intersection",
-    clipFeatures: landFC.features,
-  };
+  // Keep portion of sketch over land
 
-  // Execute one or more clip operations in order against feature
-  return clipToPolygonFeatures(feature, [keepLand], {
-    ensurePolygon: true,
-  });
+  let clipped: Feature<Polygon | MultiPolygon> | null = feature;
+
+  if (landFeatures.length === 0) {
+    clipped = null; // No land to clip to, intersection is empty
+  }
+
+  if (clipped !== null) {
+    clipped = clipMultiMerge(
+      clipped,
+      featureCollection(landFeatures),
+      "intersection",
+    );
+  }
+
+  if (!clipped || area(clipped) === 0) {
+    throw new ValidationError("Feature is outside of land boundary");
+  }
+
+  // Assume user wants the largest polygon if multiple remain
+  return biggestPolygon(clipped);
 }
 
 export default new PreprocessingHandler(clipToLand, {

packages/template-blank-project/src/functions/clipToOcean.test.ts

@@ -1,6 +1,6 @@
 import { describe, test, expect, vi, afterEach, assert } from "vitest";
 import { clipToOcean } from "./clipToOcean.js";
-import { area, featureCollection, polygon } from "@turf/turf";
+import { area, polygon } from "@turf/turf";
 
 const landFeature = polygon([
   [
@@ -54,23 +54,20 @@ describe("clipToOcean", () => {
     vi.restoreAllMocks();
   });
 
-  // Mock VectorDataSource fetchUnion method to return clipFeature
+  // Mock loadFgb method to return landFeature
   // @ts-ignore
   vi.mock(import("@seasketch/geoprocessing"), async (importOriginal) => {
     const actual = await importOriginal();
-    const VectorDataSource = vi.fn();
-    VectorDataSource.prototype.fetchUnion = vi.fn(() =>
-      featureCollection([landFeature]),
-    );
-    return { ...actual, VectorDataSource };
+    const loadFgb = vi.fn(() => [landFeature]);
+    return { ...actual, loadFgb };
   });
 
   test("clipToOcean - feature inside of land feature should throw", async () => {
     try {
       await clipToOcean(insideLand);
     } catch (error: unknown) {
       if (error instanceof Error) {
-        expect(error.message).toBe("Feature is outside of boundary");
+        expect(error.message).toBe("Feature is not in the ocean");
         return;
       }
     }

packages/template-blank-project/src/functions/clipToOcean.ts

@@ -2,18 +2,26 @@ import {
   PreprocessingHandler,
   Feature,
   Sketch,
-  clipToPolygonFeatures,
-  FeatureClipOperation,
-  VectorDataSource,
   ensureValidPolygon,
+  Polygon,
+  MultiPolygon,
+  loadFgb,
+  isPolygonFeature,
+  ValidationError,
+  clip,
+  biggestPolygon,
 } from "@seasketch/geoprocessing";
-import { bbox } from "@turf/turf";
+import { area, bbox, featureCollection } from "@turf/turf";
 
 /**
  * Preprocessor takes a Polygon feature/sketch and returns the portion that
  * is in the ocean (not on land). If results in multiple polygons then returns the largest.
  */
 export async function clipToOcean(feature: Feature | Sketch): Promise<Feature> {
+  if (!isPolygonFeature(feature)) {
+    throw new ValidationError("Input must be a polygon");
+  }
+
   // throws if not valid with specific message
   ensureValidPolygon(feature, {
     minSize: 1,
@@ -22,23 +30,27 @@ export async function clipToOcean(feature: Feature | Sketch): Promise<Feature> {
     enforceMaxSize: false,
   });
 
-  // Get land polygons - osm land vector datasource
-  const landDatasource = new VectorDataSource(
-    "https://d3p1dsef9f0gjr.cloudfront.net/",
-  );
   const featureBox = bbox(feature);
-  // one gid assigned per country, use to union subdivided pieces back together on fetch, prevents slivers
-  const landFC = await landDatasource.fetchUnion(featureBox, "gid");
 
-  const eraseLand: FeatureClipOperation = {
-    operation: "difference",
-    clipFeatures: landFC.features,
-  };
+  // Get land polygons - daylight osm land vector datasource
+  const landFeatures: Feature<Polygon | MultiPolygon>[] = await loadFgb(
+    "https://gp-global-datasources-datasets.s3.us-west-1.amazonaws.com/global-coastline-daylight-v158.fgb",
+    featureBox,
+  );
 
-  // Execute one or more clip operations in order against feature
-  return clipToPolygonFeatures(feature, [eraseLand], {
-    ensurePolygon: true,
-  });
+  // Erase portion of sketch over land
+
+  let clipped: Feature<Polygon | MultiPolygon> | null = feature;
+  if (clipped !== null && landFeatures.length > 0) {
+    clipped = clip(featureCollection([clipped, ...landFeatures]), "difference");
+  }
+
+  if (!clipped || area(clipped) === 0) {
+    throw new ValidationError("Feature is not in the ocean");
+  }
+
+  // Assume user wants the largest polygon if multiple remain
+  return biggestPolygon(clipped);
 }
 
 export default new PreprocessingHandler(clipToOcean, {

packages/template-blank-project/src/functions/clipToOceanEez.test.ts

@@ -1,6 +1,6 @@
 import { describe, test, expect, vi, afterEach, assert } from "vitest";
 import { clipToOceanEez } from "./clipToOceanEez.js";
-import { area, featureCollection, polygon } from "@turf/turf";
+import { area, polygon } from "@turf/turf";
 import { BBox } from "@seasketch/geoprocessing";
 
 const landFeature = polygon([
@@ -67,32 +67,27 @@ describe("clipToOceanEez", () => {
     vi.restoreAllMocks();
   });
 
-  // Mock VectorDataSource fetchUnion method to return clipFeature
+  // Mock loadFgb method to return landFeature
   // @ts-ignore
   vi.mock(import("@seasketch/geoprocessing"), async (importOriginal) => {
     const actual = await importOriginal();
-    const VectorDataSource = vi.fn();
-    VectorDataSource.prototype.fetchUnion = vi.fn(
-      (bbox: BBox, unionProperty?: string) => {
-        if (unionProperty === "gid") {
-          return featureCollection([landFeature]);
-        } else if (unionProperty === "UNION") {
-          return featureCollection([eezFeature]);
-        }
-      },
-    );
+    const loadFgb = vi.fn((url: string, bbox: BBox) => {
+      if (url.includes("coastline")) {
+        return [landFeature];
+      } else if (url.includes("eez")) {
+        return [eezFeature];
+      }
+    });
 
-    VectorDataSource.prototype.fetchBundle = vi.fn();
-    VectorDataSource.prototype.fetchBundleIndex = vi.fn();
-    return { ...actual, VectorDataSource };
+    return { ...actual, loadFgb };
   });
 
   test("clipToOceanEez - feature inside of land feature should throw", async () => {
     try {
       await clipToOceanEez(insideLand);
     } catch (error: unknown) {
       if (error instanceof Error) {
-        expect(error.message).toBe("Feature is outside of boundary");
+        expect(error.message).toBe("Feature is outside of EEZ boundary");
         return;
       }
     }
@@ -104,7 +99,7 @@ describe("clipToOceanEez", () => {
       await clipToOceanEez(outsideEez);
     } catch (error: unknown) {
       if (error instanceof Error) {
-        expect(error.message).toBe("Feature is outside of boundary");
+        expect(error.message).toBe("Feature is outside of EEZ boundary");
         return;
       }
     }