Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

potential fix for resamplemethod in reproject #3542

Merged
merged 3 commits into from
Oct 6, 2024
Merged

potential fix for resamplemethod in reproject #3542

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
3612b19
potential fix for resamplemethod in reproject
jdries Jun 5, 2024
dc0617b
#3541 add unit test and changelog entry
jdries Jun 10, 2024
558a654
do not catch valid error, but do avoid it in case we don't need targe…
jdries Sep 11, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
- Bump GeoTools version up to 30.x [#3521](https://github.com/locationtech/geotrellis/pull/3521)
- Dependecies update & fix "not found: type Serializable" compiler bug [#3535](https://github.com/locationtech/geotrellis/pull/3535)
- Update GDAL up to 3.9.x [#3540](https://github.com/locationtech/geotrellis/pull/3540)
- Fix reprojection with downsampling for GeotiffRasterSource and tile RDDs. [#3541](https://github.com/locationtech/geotrellis/issues/3541)

## [3.7.1] - 2024-01-08

Expand Down
20 changes: 18 additions & 2 deletions raster/src/main/scala/geotrellis/raster/reproject/RasterRegionReproject.scala
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -191,7 +191,15 @@ object RasterRegionReproject {
errorThreshold: Double
): Unit = {
val trans = Proj4Transform(dest, src)
val resampler = Resample(resampleMethod, raster.tile, raster.extent, CellSize(raster.rasterExtent.cellwidth, raster.rasterExtent.cellheight))

val targetCellSizeInSrcCRS =
try{
rasterExtent.reproject(dest,src).cellSize
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wonder if we have a broken function.

}catch {
case e:Exception => //reprojection errors happen when going outside valid area of projection, need a better fix here
raster.rasterExtent.cellSize
}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

are there any tests that cover this case?

Copy link
Contributor Author

@jdries jdries Jun 19, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

@jdries jdries Jun 19, 2024 •
edited by pomadchin
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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
    }

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hm ok I see what you need it for.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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!

Copy link
Member

@pomadchin pomadchin Sep 11, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

👍 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 👍

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

val resampler = Resample(resampleMethod, raster.tile, raster.extent, targetCellSizeInSrcCRS)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

val rowcoords = rowCoords(region, rasterExtent, trans, errorThreshold)

if (raster.cellType.isFloatingPoint) {
Expand Down Expand Up @@ -263,9 +271,17 @@ object RasterRegionReproject {
val rowcoords = rowCoords(region, rasterExtent, trans, errorThreshold)
val resampler: Array[Resample] = Array.ofDim[Resample](raster.tile.bandCount)

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
}

for (b <- 0 until raster.tile.bandCount) {
val band: Tile = raster.tile.bands(b)
resampler(b) = Resample(resampleMethod, band, raster.extent, raster.rasterExtent.cellSize)
resampler(b) = Resample(resampleMethod, band, raster.extent, targetCellSizeInSrcCRS)
}

if (raster.cellType.isFloatingPoint) {
Expand Down
20 changes: 18 additions & 2 deletions raster/src/test/scala/geotrellis/raster/reproject/ReprojectSpec.scala
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,12 +21,13 @@ import geotrellis.raster.resample._
import geotrellis.vector._
import geotrellis.raster.testkit._
import geotrellis.proj4._
import geotrellis.raster.geotiff.GeoTiffRasterSource
import geotrellis.raster.io.geotiff._

import spire.syntax.cfor._

import org.scalatest.funspec.AnyFunSpec

import java.io.File

class ReprojectSpec extends AnyFunSpec
with TileBuilders
with GeoTiffTestUtils
Expand Down Expand Up @@ -62,6 +63,21 @@ class ReprojectSpec extends AnyFunSpec
}
}

it("should (approximately) match a GDAL average interpolation on nlcd tile") {

val raster = this.createRaster(Array(1,4,1,4,1,4,1,4),4,2,CellSize(10,10))
val tempTiff = File.createTempFile("toReproject",".tif")
GeoTiff(raster, CRS.fromEpsgCode(32631)).write(tempTiff.getPath)

val rs = GeoTiffRasterSource(tempTiff.getPath).reproject(CRS.fromEpsgCode(3035), TargetCellSize(CellSize(20.0,20.0)), Average)
val reprojected = raster.reproject(CRS.fromEpsgCode(32631),CRS.fromEpsgCode(3035), Options(method = Average, errorThreshold = 0.0,targetCellSize = Some(CellSize(20,20))))

val refTile = this.createTile(Array(2.5,2.5),2,1)

assertEqual(refTile,reprojected.tile,0.1)
assertEqual(rs.read().get.tile.band(0),reprojected.tile,0.1)
}

it("should (approximately) match a GDAL nearest neighbor interpolation on slope tif") {
val raster =
SinglebandGeoTiff(geoTiffPath("reproject/slope_webmercator.tif")).raster
Expand Down