Merge pull request #28 from DanielJDufour/trimming

trim bbox reprojection, enabling global reprojection over projection bounds

Author: DanielJDufour

geowarp.js

@@ -549,7 +549,7 @@ const geowarp = function geowarp({
   let forward_turbocharged, inverse_turbocharged;
   if (turbo) {
     if (forward) {
-      out_bbox_in_srs ??= reprojectBoundingBox(out_bbox, inverse, { density: 100 });
+      out_bbox_in_srs ??= reprojectBoundingBox(out_bbox, inverse, { density: 100, nan_strategy: "skip" });
       intersect_bbox_in_srs ??= intersect(in_bbox, out_bbox_in_srs);
       forward_turbocharged = turbocharge({
         bbox: intersect_bbox_in_srs,
@@ -583,8 +583,9 @@ const geowarp = function geowarp({
   if (method === "near-vectorize" || method === "nearest-vectorize") {
     if (debug_level >= 2) console.log('[geowarp] choosing between "near" and "vectorize" for best speed');
 
-    out_bbox_in_srs ??= same_srs ? out_bbox : reprojectBoundingBox(out_bbox, inverse, { density: 100 });
+    out_bbox_in_srs ??= same_srs ? out_bbox : reprojectBoundingBox(out_bbox, inverse, { density: 100, nan_strategy: "skip" });
 
+    // average of how large each output pixel is in the input spatial reference system
     out_sample_height_in_srs = (out_bbox_in_srs[3] - out_bbox_in_srs[1]) / out_height_in_samples;
     out_sample_width_in_srs = (out_bbox_in_srs[2] - out_bbox_in_srs[0]) / out_width_in_samples;
 
@@ -613,7 +614,8 @@ const geowarp = function geowarp({
     // const [cfwd, clear_forward_cache] = cacheFunction(fwd);
 
     // reproject bounding box of output (e.g. a tile) into the spatial reference system of the input data
-    out_bbox_in_srs ??= same_srs ? out_bbox : reprojectBoundingBox(out_bbox, inverse, { density: 100 });
+    // setting nan_strategy to skip trims the box in case the output bbox extends over the bounds of the input projection
+    out_bbox_in_srs ??= same_srs ? out_bbox : reprojectBoundingBox(out_bbox, inverse, { density: 100, nan_strategy: "skip" });
     let [left, bottom, right, top] = out_bbox_in_srs;
 
     out_sample_height_in_srs ??= (top - bottom) / out_height_in_samples;
@@ -888,11 +890,18 @@ const geowarp = function geowarp({
           // combing srs reprojection and srs-to-image mapping, ensures that bounding box corners
           // are reprojected fully before calculating containing bbox
           // (prevents drift in increasing bbox twice if image is warped)
-          const [leftInRasterPixels, topInRasterPixels, rightInRasterPixels, bottomInRasterPixels] = reprojectBoundingBox(
-            [left, bottom, right, top],
-            out_srs_pt_to_in_img_pt
-          );
-
+          let leftInRasterPixels, topInRasterPixels, rightInRasterPixels, bottomInRasterPixels;
+          try {
+            [leftInRasterPixels, topInRasterPixels, rightInRasterPixels, bottomInRasterPixels] = reprojectBoundingBox(
+              [left, bottom, right, top],
+              out_srs_pt_to_in_img_pt,
+              { nan_strategy: "throw" }
+            );
+          } catch (error) {
+            // if only one pixel (or row of pixels) extends over the edge of the projection's bounds, we probably don't want to fail the whole thing
+            // an example would be warping the globe from 3857 to 4326
+            continue;
+          }
           if (debug_level >= 4) console.log("[geowarp] leftInRasterPixels:", leftInRasterPixels);
           if (debug_level >= 4) console.log("[geowarp] rightInRasterPixels:", rightInRasterPixels);
           if (debug_level >= 4) console.log("[geowarp] topInRasterPixels:", topInRasterPixels);

test.js

@@ -860,3 +860,67 @@ test("antarctica with NaN", async ({ eq }) => {
     }
   }
 });
+
+test("issue 27: globe 3857 to 4326", async ({ eq }) => {
+  const filename = "gadas-world.tif";
+  const filepath = path.resolve(__dirname, "./test-data", filename);
+  const geotiff = await GeoTIFF.fromFile(filepath);
+  const image = await geotiff.getImage(0);
+  const in_data = await image.readRasters();
+  const bbox = getBoundingBox(image);
+  const in_height = image.getHeight();
+  const in_width = image.getWidth();
+  const out_height = 180 * 2;
+  const out_width = 360 * 2;
+
+  const { forward, inverse } = proj4("EPSG:3857", "EPSG:4326");
+
+  const methods = ["near-vectorize", "near", "bilinear", "median"];
+  const turbos = [true, false];
+  for (let i = 0; i < methods.length; i++) {
+    const method = methods[i];
+    for (let ii = 0; ii < 2; ii++) {
+      const turbo = turbos[ii];
+      const result = await geowarp({
+        debug_level: 0,
+        in_bbox: bbox,
+        in_geotransform: [-20057076.25595305, 39135.75848200009, 0, 12640208.839021027, 0, -39135.75848200009],
+        in_data,
+        in_layout: "[band][row,column]",
+        in_srs: 3857,
+        in_height,
+        in_width,
+        out_bbox: [-180, -90, 180, 90],
+        out_height,
+        out_no_data: null,
+        out_width,
+        out_layout: "[band][row][column]",
+        out_srs: 4326,
+        forward,
+        inverse,
+        method,
+        round: true,
+        theoretical_max: 255,
+        theoretical_min: 0,
+        turbo
+      });
+      // console.log(method + " warped");
+      // console.dir(result.data, { maxArrayLength: 3 });
+
+      // console.dir(
+      //   result.data.map(band => band.map(row => row[0])),
+      //   { maxArrayLength: 500 }
+      // );
+
+      // check that no NaN values in output
+      eq(
+        result.data.flat(3).findIndex(it => isNaN(it)),
+        -1
+      );
+
+      if (process.env.WRITE) {
+        writePNGSync({ h: out_height, w: out_width, data: result.data, filepath: `./test-output/gadas-whole-4326-${method}${turbo ? "-turbo" : ""}` });
+      }
+    }
+  }
+});

test.ts

@@ -377,6 +377,7 @@ const runTileTests = async ({
       "25,33,20",
       "27,35,22",
       "28,30,17",
+      "32,34,21",
       "33,43,34",
       "36,46,45",
       "40,49,47",