Fix MosaicRasterSource, GDALRasterSource and GeoTiffResampleRasterSource behavior

CHANGELOG.md

@@ -20,6 +20,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - GeoTrellisRasterSource should return None on empty reads [#3240](https://github.com/locationtech/geotrellis/pull/3240) 
 - ArrayTile equals method always returns true if first elements are NaN [#3242](https://github.com/locationtech/geotrellis/issues/3242)
 - Fixed resource issue with JpegDecompressor that was causing a "too many open files in the system" exception on many parallel reads of JPEG compressed GeoTiffs. [#3249](https://github.com/locationtech/geotrellis/pull/3249)
+- Fix MosaicRasterSource, GDALRasterSource, GeoTiffResample raster sources behavior  [#3252](https://github.com/locationtech/geotrellis/pull/3252)
 
 ## [3.3.0] - 2020-04-07
 

gdal/src/main/scala/geotrellis/raster/gdal/GDALWarpOptions.scala

@@ -256,14 +256,32 @@ case class GDALWarpOptions(
    */
   def reproject(rasterExtent: GridExtent[Long], sourceCRS: CRS, targetCRS: CRS, resampleTarget: ResampleTarget = DefaultTarget, resampleMethod: ResampleMethod = NearestNeighbor): GDALWarpOptions = {
     val reprojectOptions = ResampleTarget.toReprojectOptions(rasterExtent, resampleTarget, resampleMethod)
-    val re = rasterExtent.reproject(sourceCRS, targetCRS, reprojectOptions)
+    val reprojectedRasterExtent = rasterExtent.reproject(sourceCRS, targetCRS, reprojectOptions)
 
-    this.copy(
-      cellSize       = re.cellSize.some,
-      targetCRS      = targetCRS.some,
-      sourceCRS      = sourceCRS.some,
-      resampleMethod = reprojectOptions.method.some
-    )
+    resampleTarget match {
+      case TargetDimensions(cols, rows) =>
+        this.copy(
+          te             = reprojectedRasterExtent.extent.some,
+          cellSize       = None,
+          targetCRS      = targetCRS.some,
+          sourceCRS      = this.sourceCRS orElse sourceCRS.some,
+          resampleMethod = reprojectOptions.method.some,
+          dimensions     = (cols.toInt, rows.toInt).some
+        )
+      case _ =>
+        val re = {
+          val targetRasterExtent = resampleTarget(reprojectedRasterExtent).toRasterExtent
+          if(this.alignTargetPixels) targetRasterExtent.alignTargetPixels else targetRasterExtent
+        }
+
+        this.copy(
+          cellSize       = re.cellSize.some,
+          te             = re.extent.some,
+          targetCRS      = targetCRS.some,
+          sourceCRS      = this.sourceCRS orElse sourceCRS.some,
+          resampleMethod = reprojectOptions.method.some
+        )
+    }
   }
 
   /** Adjust GDAL options to represents resampling with following parameters .

gdal/src/test/scala/geotrellis/raster/gdal/GDALRasterSourceSpec.scala

@@ -16,15 +16,15 @@
 
 package geotrellis.raster.gdal
 
-import geotrellis.proj4.LatLng
+import cats.data.NonEmptyList
+import geotrellis.proj4.{CRS, LatLng, WebMercator}
 import geotrellis.raster.geotiff.GeoTiffRasterSource
 import geotrellis.raster._
-import geotrellis.raster.io.geotiff.GeoTiff
+import geotrellis.raster.io.geotiff.{AutoHigherResolution, GeoTiff}
 import geotrellis.raster.io.geotiff.reader.GeoTiffReader
 import geotrellis.raster.resample._
 import geotrellis.raster.testkit._
 import geotrellis.vector.Extent
-
 import org.scalatest._
 
 class GDALRasterSourceSpec extends FunSpec with RasterMatchers with GivenWhenThen {
@@ -114,6 +114,45 @@ class GDALRasterSourceSpec extends FunSpec with RasterMatchers with GivenWhenThe
       source.metadata.attributes.mapValues(_.toUpperCase) shouldBe tsource.metadata.attributes.mapValues(_.toUpperCase)
     }
 
+    it("should perform a chained reprojection") {
+      val rs = GDALRasterSource(Resource.path("vlm/lc8-utm-1.tif"))
+      val drs = rs.reproject(WebMercator).reproject(LatLng).asInstanceOf[GDALRasterSource]
+      val rrs = rs.reproject(LatLng).asInstanceOf[GDALRasterSource]
+
+      drs.crs shouldBe rrs.crs
+
+      val RasterExtent(Extent(dxmin, dymin, dxmax, dymax), dcw, dch, dc, dr) = drs.gridExtent.toRasterExtent()
+      val RasterExtent(Extent(rxmin, rymin, rxmax, rymax), rcw, rch, rc, rr) = rrs.gridExtent.toRasterExtent()
+
+      // our reprojection rounding error
+      // to fix it we can round cellSize in the options reprojection by 6 numbers after the decimal point
+      dxmin shouldBe rxmin +- 1e-1
+      dymin shouldBe rymin +- 1e-1
+      dxmax shouldBe rxmax +- 1e-1
+      dymax shouldBe rymax +- 1e-1
+
+      dcw shouldBe rcw +- 1e-6
+      dch shouldBe rch +- 1e-6
+
+      dc shouldBe rc +- 1
+      dr shouldBe rr +- 1
+
+      // manually ensure that the difference in two GDALWarpOptions
+      drs.options.copy(cellSize = None, te = None) shouldBe rrs.options.copy(cellSize = None, te = None)
+      val CellSize(docw, doch) = drs.options.cellSize.get
+      val CellSize(rocw, roch) = rrs.options.cellSize.get
+      docw shouldBe rocw +- 1e-6
+      doch shouldBe roch +- 1e-6
+
+      val Extent(dtexmin, dteymin, dtexmax, dteymax) = drs.options.te.get
+      val Extent(rtexmin, rteymin, rtexmax, rteymax) = rrs.options.te.get
+
+      dtexmin shouldBe rtexmin +- 1e-1
+      dteymin shouldBe rteymin +- 1e-1
+      dtexmax shouldBe rtexmax +- 1e-1
+      dteymax shouldBe rteymax +- 1e-1
+    }
+
     describe("should derive the cellType consistently with GeoTiffRasterSource") {
       val raster = Raster(
         ArrayTile(Array(
@@ -158,5 +197,39 @@ class GDALRasterSourceSpec extends FunSpec with RasterMatchers with GivenWhenThe
         }
       }
     }
+
+
+    it("reprojectToRegion should produce an expected raster") {
+      val expected = GeoTiffRasterSource(Resource.path("vlm/lc8-utm-re-mosaic-expected.tif")).read().get
+
+      val uris = List(
+        Resource.path("vlm/lc8-utm-1.tif"),
+        Resource.path("vlm/lc8-utm-2.tif"),
+        Resource.path("vlm/lc8-utm-3.tif")
+      )
+
+      val targetGridExtent = GridExtent[Long](
+        Extent(-89.18876450968908, 37.796833507913995, -73.9486554460876, 41.41061537114088),
+        CellSize(0.01272458402544677,0.01272129304140357)
+      )
+      val e = Extent(-89.18876450968908, 37.796833507913995, -73.9486554460876, 41.41061537114088)
+      val targetCRS = CRS.fromEpsgCode(4326)
+
+      val rasterSources = NonEmptyList.fromListUnsafe(uris.map(GDALRasterSource(_)))
+      val moisac = MosaicRasterSource.instance(rasterSources, rasterSources.head.crs)
+
+      val actual =
+        moisac
+          .reprojectToRegion(
+            targetCRS,
+            targetGridExtent.toRasterExtent,
+            NearestNeighbor,
+            AutoHigherResolution
+          )
+          .read(e)
+          .get
+
+      assertEqual(actual, expected)
+    }
   }
 }

raster/src/main/scala/geotrellis/raster/MosaicRasterSource.scala

@@ -19,7 +19,6 @@ package geotrellis.raster
 import geotrellis.vector._
 import geotrellis.raster._
 import geotrellis.raster.resample._
-import geotrellis.raster.reproject.Reproject
 import geotrellis.proj4.{CRS, WebMercator}
 import geotrellis.raster.io.geotiff.OverviewStrategy
 
@@ -86,13 +85,8 @@ abstract class MosaicRasterSource extends RasterSource {
     resampleTarget: ResampleTarget = DefaultTarget,
     method: ResampleMethod = ResampleMethod.DEFAULT,
     strategy: OverviewStrategy = OverviewStrategy.DEFAULT
-  ): RasterSource =
-    MosaicRasterSource(
-      sources map { _.reproject(targetCRS, resampleTarget, method, strategy) },
-      crs,
-      gridExtent.reproject(this.crs, targetCRS, Reproject.Options.DEFAULT.copy(method = method)),
-      name
-    )
+   ): RasterSource =
+    MosaicRasterSource.instance(sources.map(_.reproject(targetCRS, resampleTarget, method, strategy)), targetCRS, name)
 
   def read(extent: Extent, bands: Seq[Int]): Option[Raster[MultibandTile]] = {
     val rasters = sources map { _.read(extent, bands) }
@@ -108,11 +102,13 @@ abstract class MosaicRasterSource extends RasterSource {
     resampleTarget: ResampleTarget,
     method: ResampleMethod,
     strategy: OverviewStrategy
-  ): RasterSource = MosaicRasterSource(
-    sources map { _.resample(resampleTarget, method, strategy) }, crs, name)
+  ): RasterSource =
+    MosaicRasterSource.instance(sources.map(_.resample(resampleTarget, method, strategy)), crs, name)
 
   def convert(targetCellType: TargetCellType): RasterSource =
-    MosaicRasterSource(sources map { _.convert(targetCellType) }, crs, name)
+    MosaicRasterSource.instance(sources.map(_.convert(targetCellType)), crs, name)
+
+  override def toString: String = s"MosaicRasterSource(${sources.toList}, $crs, $gridExtent, $name)"
 }
 
 object MosaicRasterSource {
@@ -123,7 +119,7 @@ object MosaicRasterSource {
       def combine(l: Raster[MultibandTile], r: Raster[MultibandTile]) = {
         val targetRE = RasterExtent(
           l.rasterExtent.extent combine r.rasterExtent.extent,
-          List(l.rasterExtent.cellSize, r.rasterExtent.cellSize).minBy(_.resolution))
+          List(l.rasterExtent.cellSize, r.rasterExtent.cellSize).maxBy(_.resolution))
         val result = l.resample(targetRE) merge r.resample(targetRE)
         result
       }
@@ -144,13 +140,34 @@ object MosaicRasterSource {
       }
     }
 
+  def instance(sourcesList: NonEmptyList[RasterSource], targetCRS: CRS, targetGridExtent: GridExtent[Long], sourceName: SourceName): MosaicRasterSource =
+    new MosaicRasterSource {
+      val sources: NonEmptyList[RasterSource] = sourcesList
+      val crs: CRS = targetCRS
+      def gridExtent: GridExtent[Long] = targetGridExtent
+      val name: SourceName = sourceName
+    }
+
+  def instance(sourcesList: NonEmptyList[RasterSource], targetCRS: CRS, sourceName: SourceName): MosaicRasterSource = {
+    val combinedExtent = sourcesList.map(_.extent).toList.reduce(_ combine _)
+    val minCellSize = sourcesList.map(_.cellSize).toList.maxBy(_.resolution)
+    val combinedGridExtent = GridExtent[Long](combinedExtent, minCellSize)
+    instance(sourcesList, targetCRS, combinedGridExtent, sourceName)
+  }
+
+  def instance(sourcesList: NonEmptyList[RasterSource], targetCRS: CRS): MosaicRasterSource =
+    instance(sourcesList, targetCRS, EmptyName)
+
+  def instance(sourcesList: NonEmptyList[RasterSource], targetCRS: CRS, targetGridExtent: GridExtent[Long]): MosaicRasterSource =
+    instance(sourcesList, targetCRS, targetGridExtent, EmptyName)
+
   def apply(sourcesList: NonEmptyList[RasterSource], targetCRS: CRS, targetGridExtent: GridExtent[Long]): MosaicRasterSource =
     apply(sourcesList, targetCRS, targetGridExtent, EmptyName)
 
   def apply(sourcesList: NonEmptyList[RasterSource], targetCRS: CRS, targetGridExtent: GridExtent[Long], rasterSourceName: SourceName): MosaicRasterSource = {
     new MosaicRasterSource {
       val name = rasterSourceName
-      val sources = sourcesList map { _.reprojectToGrid(targetCRS, gridExtent) }
+      val sources = sourcesList.map(_.reprojectToGrid(targetCRS, gridExtent))
       val crs = targetCRS
 
       def gridExtent: GridExtent[Long] = targetGridExtent
@@ -166,16 +183,10 @@ object MosaicRasterSource {
       val sources = sourcesList map { _.reprojectToGrid(targetCRS, sourcesList.head.gridExtent) }
       val crs = targetCRS
       def gridExtent: GridExtent[Long] = {
-        val reprojectedExtents =
-          sourcesList map { source => source.gridExtent.reproject(source.crs, targetCRS) }
-        val minCellSize: CellSize = reprojectedExtents.toList map { rasterExtent =>
-          CellSize(rasterExtent.cellwidth, rasterExtent.cellheight)
-        } minBy { _.resolution }
-        reprojectedExtents.toList.reduce(
-          (re1: GridExtent[Long], re2: GridExtent[Long]) => {
-            re1.withResolution(minCellSize) combine re2.withResolution(minCellSize)
-          }
-        )
+        val reprojectedSources = sources.toList
+        val combinedExtent = reprojectedSources.map(_.extent).reduce(_ combine _)
+        val minCellSize = reprojectedSources.map(_.cellSize).maxBy(_.resolution)
+        GridExtent[Long](combinedExtent, minCellSize)
       }
     }
   }

raster/src/main/scala/geotrellis/raster/RasterSource.scala

@@ -73,7 +73,7 @@ abstract class RasterSource extends CellGrid[Long] with RasterMetadata {
     */
   def reprojectToRegion(targetCRS: CRS, region: RasterExtent, method: ResampleMethod = ResampleMethod.DEFAULT, strategy: OverviewStrategy = OverviewStrategy.DEFAULT): RasterSource =
     if (targetCRS == this.crs && region == this.gridExtent) this
-    else if (targetCRS == this.crs) resampleToRegion(region.asInstanceOf[GridExtent[Long]], method, strategy)
+    else if (targetCRS == this.crs) resampleToRegion(region.toGridType[Long], method, strategy)
     else reprojection(targetCRS, TargetRegion(region.toGridType[Long]), method, strategy)
 
   def resample(resampleTarget: ResampleTarget, method: ResampleMethod, strategy: OverviewStrategy): RasterSource

raster/src/main/scala/geotrellis/raster/geotiff/GeoTiffReprojectRasterSource.scala

@@ -71,13 +71,7 @@ class GeoTiffReprojectRasterSource(
         Reproject.Options.DEFAULT.copy(method = resampleMethod, errorThreshold = errorThreshold)
       )
 
-    resampleTarget match {
-      case targetRegion: TargetRegion => targetRegion.region.toGridType[Long]
-      case targetAlignment: TargetAlignment => targetAlignment(reprojectedRasterExtent)
-      case targetDimensions: TargetDimensions => targetDimensions(reprojectedRasterExtent)
-      case targetCellSize: TargetCellSize => targetCellSize(reprojectedRasterExtent)
-      case _ => reprojectedRasterExtent
-    }
+    resampleTarget(reprojectedRasterExtent)
   }
 
   lazy val resolutions: List[CellSize] =

raster/src/main/scala/geotrellis/raster/geotiff/GeoTiffResampleRasterSource.scala

@@ -57,6 +57,7 @@ class GeoTiffResampleRasterSource(
 
   def crs: CRS = tiff.crs
 
+  // GDAL aligns target pixels, do we need to do it as well?
   override lazy val gridExtent: GridExtent[Long] = resampleTarget(tiff.rasterExtent.toGridType[Long])
 
   lazy val resolutions: List[CellSize] = tiff.cellSize :: tiff.overviews.map(_.cellSize)
@@ -74,13 +75,7 @@ class GeoTiffResampleRasterSource(
             Reproject.Options.DEFAULT.copy(method = resampleMethod, errorThreshold = errorThreshold)
           )
 
-        resampleTarget match {
-          case targetRegion: TargetRegion => targetRegion.region.toGridType[Long]
-          case targetAlignment: TargetAlignment => targetAlignment(reprojectedRasterExtent)
-          case targetDimensions: TargetDimensions => targetDimensions(reprojectedRasterExtent)
-          case targetCellSize: TargetCellSize => targetCellSize(reprojectedRasterExtent)
-          case _ => reprojectedRasterExtent
-        }
+        resampleTarget(reprojectedRasterExtent)
       }
     }
 
@@ -124,10 +119,7 @@ class GeoTiffResampleRasterSource(
 
     geoTiffTile.crop(windows.keys.toSeq, bands.toArray).map { case (gb, tile) =>
       val targetRasterExtent = windows(gb)
-      Raster(
-        tile = tile,
-        extent = targetRasterExtent.extent
-      ).resample(targetRasterExtent.cols, targetRasterExtent.rows, method)
+      Raster(tile, closestTiffOverview.rasterExtent.extentFor(gb, clamp = true)).resample(targetRasterExtent, method)
     }
   }