potential fix for resamplemethod in reproject (#3542)

* address resamplemethod in reproject #3541
locationtech · Oct 6, 2024 · cfa8e76 · cfa8e76
1 parent 3b53c62
commit cfa8e76
Show file tree

Hide file tree

Showing 3 changed files with 32 additions and 4 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -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. Reprojection outside the valid projection bounds may now throw a GeoAttrsError. [#3541](https://github.com/locationtech/geotrellis/issues/3541)
 
 ## [3.7.1] - 2024-01-08
 

diff --git a/raster/src/main/scala/geotrellis/raster/reproject/RasterRegionReproject.scala b/raster/src/main/scala/geotrellis/raster/reproject/RasterRegionReproject.scala
@@ -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 resampler = resampleMethod match {
+        case resampleMethod1: PointResampleMethod =>
+          Resample(resampleMethod1, raster.tile, raster.extent)
+        case _ =>
+          //throws GeoAttrsError when applied for invalid extent
+          val targetCellSizeInSrcCRS = rasterExtent.reproject(dest, src).cellSize
+          Resample(resampleMethod, raster.tile, raster.extent, targetCellSizeInSrcCRS)
+      }
       val rowcoords = rowCoords(region, rasterExtent, trans, errorThreshold)
 
       if (raster.cellType.isFloatingPoint) {
@@ -263,9 +271,12 @@ object RasterRegionReproject {
       val rowcoords = rowCoords(region, rasterExtent, trans, errorThreshold)
       val resampler: Array[Resample] = Array.ofDim[Resample](raster.tile.bandCount)
 
+      val targetCellSizeInSrcCRS = rasterExtent.reproject(dest,src).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) {

diff --git a/raster/src/test/scala/geotrellis/raster/reproject/ReprojectSpec.scala b/raster/src/test/scala/geotrellis/raster/reproject/ReprojectSpec.scala
@@ -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
@@ -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