potential fix for resamplemethod in reproject #3542

jdries · 2024-06-05T08:30:40Z

Overview

Proposed fix for a very important bug in reprojecting to lower resolution.
I think it's better (more correct) than previous behaviour, but not sure yet if it is 100% accurate.

It seems to me that the downsampling now happens in the original projection system, and then it is warping to the new CRS. This will work fine if the low resolution pixel is rectangular in both projection systems.

Checklist

./CHANGELOG.md updated, if necessary. Link to the issue if closed, otherwise the PR.
Module Hierarchy updated, if necessary
docs guides update, if necessary
New user API has useful Scaladoc strings
Unit tests added for bug-fix or new feature

Demo

Optional. Screenshots/REPL

Notes

Optional. Ancillary topics, caveats, alternative strategies that didn't work out, anything else.

Closes #XXX

locationtech#3541

jdries · 2024-06-10T12:54:59Z

The fix works but I had to fall back to old (bad) behaviour in cases at the edge of valid target crs bounds. My method to compute target cellsize resulted in an error there.

pomadchin · 2024-06-15T11:17:23Z

raster/src/main/scala/geotrellis/raster/reproject/RasterRegionReproject.scala

+      val targetCellSizeInSrcCRS =
+        try{
+          rasterExtent.reproject(dest,src).cellSize
+        }catch {
+          case e:Exception => //reprojection errors happen when going outside valid area of projection, need a better fix here
+            raster.rasterExtent.cellSize
+        }


are there any tests that cover this case?

Not yet, but I now had some time to dig deeper.
The original stack trace was:

at geotrellis.raster.GridExtent.<init>(GridExtent.scala:46) at geotrellis.raster.GridExtent.<init>(GridExtent.scala:58) at geotrellis.raster.reproject.ReprojectRasterExtent$.apply(ReprojectRasterExtent.scala:71) at geotrellis.raster.GridExtent$gridExtentMethods.$anonfun$reproject$2(GridExtent.scala:432) at scala.Option.getOrElse(Option.scala:189) at geotrellis.raster.GridExtent$gridExtentMethods.reproject(GridExtent.scala:432) at geotrellis.raster.GridExtent$gridExtentMethods.reproject(GridExtent.scala:435) at org.openeo.geotrellis.reproject.RasterRegionReproject$$anon$2.reprojectToBuffer(RasterRegionReproject.scala:268) at org.openeo.geotrellis.reproject.RasterRegionReproject$$anon$2.regionReproject(RasterRegionReproject.scala:229) at org.openeo.geotrellis.reproject.TileRDDReproject$.createCombiner$1(TileRDDReproject.scala:192) at org.openeo.geotrellis.reproject.TileRDDReproject$.$anonfun$apply$6(TileRDDReproject.scala:206)

and got triggered by a more complicated test case here:
https://github.com/Open-EO/openeo-geotrellis-extensions/blob/106784e1cb83a87f1bacef3beb1b9df77c5ecb65/openeo-geotrellis/src/test/scala/org/openeo/geotrellis/OpenEOProcessesSpec.scala#L582

I could reproduce it in a simpler test, getting a similar but not exactly the same stack trace:

it("should (approximately) match a GDAL average interpolation on nlcd tile at edge of valid bounds") { val tile = this.createTile(Array(1,4,1,4,1,4,1,4),4,2) val raster = Raster(tile,Extent(-40.50000387499999, -89.99999, 40.499993875000015, -51.749994249999986)) val tempTiff = File.createTempFile("toReproject",".tif") GeoTiff(raster, LatLng).write(tempTiff.getPath) val tiffRe = RasterExtent(Extent(-4452779.631730943, -6.325724963354129E7, 0.0, -5.839130735403811E7), CellSize(15000.0,15000.0)) val alignment = TargetRegion(tiffRe) val rs = GeoTiffReprojectRasterSource(tempTiff.getPath,WebMercator, alignment, Average) val reprojectedFromRS = rs.read().get.tile.band(0) val reprojected = raster.reproject(LatLng,WebMercator, Options(method = Average, errorThreshold = 0.0,targetCellSize = Some(CellSize(1000,1000)))) val refTile = this.createTile(Array(2.5,2.5),2226, 10363) //assertEqual(refTile,reprojected.tile,0.1) assertEqual(reprojectedFromRS,reprojected.tile,0.1) }

is giving:

invalid rows: 0 geotrellis.raster.GeoAttrsError: invalid rows: 0 at geotrellis.raster.GridExtent.<init>(GridExtent.scala:46) at geotrellis.raster.GridExtent.<init>(GridExtent.scala:58) at geotrellis.raster.reproject.ReprojectRasterExtent$.apply(ReprojectRasterExtent.scala:71) at geotrellis.raster.reproject.ReprojectRasterExtent$.apply(ReprojectRasterExtent.scala:76) at geotrellis.raster.geotiff.GeoTiffReprojectRasterSource.closestTiffOverview$lzycompute(GeoTiffReprojectRasterSource.scala:86) at geotrellis.raster.geotiff.GeoTiffReprojectRasterSource.closestTiffOverview(GeoTiffReprojectRasterSource.scala:81) at geotrellis.raster.geotiff.GeoTiffReprojectRasterSource.readBounds(GeoTiffReprojectRasterSource.scala:110) at geotrellis.raster.geotiff.GeoTiffReprojectRasterSource.read(GeoTiffReprojectRasterSource.scala:98) at geotrellis.raster.geotiff.GeoTiffReprojectRasterSource.read(GeoTiffReprojectRasterSource.scala:94) at geotrellis.raster.RasterSource.read(RasterSource.scala:134) at geotrellis.raster.reproject.ReprojectSpec.$anonfun$new$5(ReprojectSpec.scala:92)

What's happening here is that we are going across valid bounds of webmercator. The root cause is then in these lines:
https://github.com/VitoTAP/geotrellis/blob/b071b333581ed7206ad14a06b1b70cfde6fb65a7/raster/src/main/scala/geotrellis/raster/reproject/ReprojectRasterExtent.scala#L63

The 'toLong' will round down to zero, which is not accepted later on in the code.
When that case happens, I think it also doesn't really matter any more what cellsize is used, because we warp into an area with cols or rows == 0.

ic ic; the thing is that this try / catch is a bit sketchy 🤔 we definitely need to find a better way addressing it / or at least make exception more specific.

Here's a case that really isolates and exactly reproduces the issue. Probably these extents are fairly problematic to begin with:

it("should (approximately) match a GDAL average interpolation on nlcd tile at edge of valid bounds") { val tile = this.createNoData(288,272,ByteConstantNoDataCellType) val raster = Raster(tile,Extent(-40.50000387499999, -89.99999, 40.499993875000015, -51.749994249999986)) val tiffRe = RasterExtent(Extent(-4452779.631730943, -6.325724963354129E7, 0.0, -5.839130735403811E7), CellSize(14744.303416327626,16112.391653984045)) val rr = implicitly[RasterRegionReproject[Tile]] val reprojectedRasterRegion: Raster[Tile] = rr.regionReproject( raster, LatLng, WebMercator, tiffRe, tiffRe.extent.toPolygon(), Average, 0.2 ) //succesfull if it doesn't throw an exception }

Hm ok I see what you need it for.

I think the change is good, sadly I don't think we can merge this workaround with throws; it simply uses incorrect resolutions as an input for the functions.

Mb it is acceptable (it worked somehow up until now) but not sure that's what we really need to do here 🤔 The error thrown seems to be legit.

Hi @pomadchin, the error is indeed legit for sure. The biggest problem is that if we start throwing it, some existing code might break, because now it is perfectly possible to do reprojection outside the valid bounds of the coordinate system.
However, from my side, this would be fine, it's at least better than having resampling fully broken. So shall I update the PR to simply not catch the exception, so that downstream users can deal with it?
Whatever your preferred way forward is, I'd very much like to fix the resampling!

👍 yes, the change is very good, that's an aweome bug catch.

The error handling via try is no good though, it is expected to fail 👍

If you could remove the error catching we can merge it in I think.

committed it, I did still came up with a clean way to keep the 'old' behaviour intact when using a resample method that does not require computing this cellsize. This will reduce upgrade issues for a lot of users.

pomadchin · 2024-06-30T13:42:06Z

raster/src/main/scala/geotrellis/raster/reproject/RasterRegionReproject.scala

+          case e:Exception => //reprojection errors happen when going outside valid area of projection, need a better fix here
+            raster.rasterExtent.cellSize
+        }
+      val resampler = Resample(resampleMethod, raster.tile, raster.extent, targetCellSizeInSrcCRS)


I see, the issue is that the targetCellSize was / is in the wrong CRS.

pomadchin · 2024-06-30T13:42:21Z

raster/src/main/scala/geotrellis/raster/reproject/RasterRegionReproject.scala

+
+      val targetCellSizeInSrcCRS =
+        try{
+          rasterExtent.reproject(dest,src).cellSize


I wonder if we have a broken function.

…t cellsize

jdries · 2024-09-19T08:20:23Z

@pomadchin ready to merge?

pomadchin · 2024-09-19T13:07:24Z

@jdries I'll take a look this week!

pomadchin · 2024-10-06T16:09:59Z

Thank you!

jdries added 2 commits June 5, 2024 10:27

potential fix for resamplemethod in reproject

3612b19

locationtech#3541

locationtech#3541 add unit test and changelog entry

dc0617b

pomadchin reviewed Jun 15, 2024

View reviewed changes

pomadchin added the bug label Jun 19, 2024

pomadchin reviewed Jun 30, 2024

View reviewed changes

jdries mentioned this pull request Sep 9, 2024

Expose area of use locationtech/proj4j#113

Open

pomadchin self-requested a review September 11, 2024 02:37

pomadchin assigned jdries Sep 11, 2024

do not catch valid error, but do avoid it in case we don't need targe…

558a654

…t cellsize

pomadchin merged commit cfa8e76 into locationtech:master Oct 6, 2024
6 checks passed

pomadchin mentioned this pull request Oct 6, 2024

RDD reproject does not correctly apply resampleMethod #3541

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

potential fix for resamplemethod in reproject #3542

potential fix for resamplemethod in reproject #3542

jdries commented Jun 5, 2024 •

edited

Loading

jdries commented Jun 10, 2024

pomadchin Jun 15, 2024

jdries Jun 19, 2024 •

edited

Loading

pomadchin Jun 19, 2024

jdries Jun 19, 2024 •

edited by pomadchin

Loading

pomadchin Jun 30, 2024

pomadchin Jun 30, 2024

jdries Sep 9, 2024

pomadchin Sep 11, 2024 •

edited

Loading

pomadchin Sep 11, 2024

jdries Sep 11, 2024

pomadchin Jun 30, 2024

pomadchin Jun 30, 2024

jdries commented Sep 19, 2024

pomadchin commented Sep 19, 2024

pomadchin commented Oct 6, 2024

potential fix for resamplemethod in reproject #3542

potential fix for resamplemethod in reproject #3542

Conversation

jdries commented Jun 5, 2024 • edited Loading

Overview

Checklist

Demo

Notes

jdries commented Jun 10, 2024

Choose a reason for hiding this comment

jdries Jun 19, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jdries Jun 19, 2024 • edited by pomadchin Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pomadchin Sep 11, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jdries commented Sep 19, 2024

pomadchin commented Sep 19, 2024

pomadchin commented Oct 6, 2024

jdries commented Jun 5, 2024 •

edited

Loading

jdries Jun 19, 2024 •

edited

Loading

jdries Jun 19, 2024 •

edited by pomadchin

Loading

pomadchin Sep 11, 2024 •

edited

Loading