From 2e52a1c3627d14eb6873eca2932fe94f6b14c470 Mon Sep 17 00:00:00 2001 From: DanielJDufour Date: Sat, 24 Feb 2024 19:17:07 -0500 Subject: [PATCH] trim bbox reprojection, enabling global reprojection over projection bounds --- geowarp.js | 25 ++++++++++++++------- test.js | 64 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ test.ts | 1 + 3 files changed, 82 insertions(+), 8 deletions(-) diff --git a/geowarp.js b/geowarp.js index 8d3d625..ae9c551 100644 --- a/geowarp.js +++ b/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); diff --git a/test.js b/test.js index 6df4d6a..4213c5f 100644 --- a/test.js +++ b/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" : ""}` }); + } + } + } +}); diff --git a/test.ts b/test.ts index a233f83..6829739 100644 --- a/test.ts +++ b/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",