Fix MosaicRasterSource.tileToLayout behavior

CHANGELOG.md

@@ -16,6 +16,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - GDALWarpOptions incorrect rpc flag [#3370](https://github.com/locationtech/geotrellis/issues/3370)
 - S3LayerDeleter cannot handle over 1000 objects to delete [#3371](https://github.com/locationtech/geotrellis/issues/3371)
 - Drop Scala 2.11 cross compilation [#3259](https://github.com/locationtech/geotrellis/issues/3259)
+- Fix MosaicRasterSource.tileToLayout behavior [#3338](https://github.com/locationtech/geotrellis/pull/3338)
 
 ## [3.5.2] - 2021-02-01
 

layer/src/main/scala/geotrellis/layer/LayoutTileSource.scala

@@ -82,10 +82,22 @@ class LayoutTileSource[K: SpatialComponent](
       if (raster.tile.cols == layout.tileCols && raster.tile.rows == layout.tileRows) {
         raster.tile
       } else {
-        // raster is smaller but not bigger than I think ...
-        // its offset is relative to the raster we wished we had
-        val colOffset = bounds.colMin - sourcePixelBounds.colMin
-        val rowOffset = bounds.rowMin - sourcePixelBounds.rowMin
+        /** Raster is smaller but not bigger than I think ...
+          * Its offset is relative to the raster we wished we had.
+          * When the source represents a single [[RasterSource]], it is possible co compute offsets this way:
+          *   val colOffset = bounds.colMin - sourcePixelBounds.colMin
+          *   val rowOffset = bounds.rowMin - sourcePixelBounds.rowMin
+          *
+          * However, if the source is a [[MosaicRasterSource]], offsets are still relative to the underlying [[RasterSource]].
+          * In this case it is possible to compute offsets through the returned extent.
+          */
+
+        // the actual returned extent
+        val actualExtent = raster.extent
+        // expected extent
+        val expectedRasterExtent = RasterExtent(layout.mapTransform(spatialComponent), layout.tileCols, layout.tileRows)
+        val (colOffset, rowOffset) = expectedRasterExtent.mapToGrid(actualExtent.xmin, actualExtent.ymax)
+
         raster.tile.mapBands { (_, band) =>
           PaddedTile(band, colOffset.toInt, rowOffset.toInt, layout.tileCols, layout.tileRows)
         }

layer/src/test/scala/geotrellis/layer/LayoutTileSourceSpec.scala

@@ -23,12 +23,15 @@ import geotrellis.raster.geotiff._
 import geotrellis.raster.io.geotiff.reader._
 import geotrellis.raster.resample._
 import geotrellis.vector._
+import geotrellis.raster.io.geotiff.AutoHigherResolution
 
-import java.io.File
+import cats.data.NonEmptyList
 
+import java.io.File
 import org.scalatest.funspec.AnyFunSpec
 
 class LayoutTileSourceSpec extends AnyFunSpec with RasterMatchers {
+
   /** Function to get paths from the raster/data dir.  */
   def rasterGeoTiffPath(name: String): String = {
     def baseDataPath = "raster/data"
@@ -268,13 +271,54 @@ class LayoutTileSourceSpec extends AnyFunSpec with RasterMatchers {
           .get
           .band(0)
         val arr = tile.toArray
-        val ones = tile.mapIfSet(i => 1).toArray()
-        ones.sum shouldBe(arr.size)
+        val ones = tile.mapIfSet(_ => 1).toArray()
+        ones.sum shouldBe arr.size
       }
 
       neighborhood map {
         case (x, y) => checkExtent(x, y, rs)
       }
     }
+
+    // https://github.com/locationtech/geotrellis/issues/3253
+    it("MosaicRasterSource should be properly tiled to layout") {
+      val ld = LayoutDefinition(
+        Extent(-2.003750834278925E7, -2.003750834278925E7, 2.003750834278925E7, 2.003750834278925E7),
+        TileLayout(256, 256, 256, 256)
+      )
+
+      // For SpatialKey(74, 96):
+      // rs1: GridBounds(133,0,388,12)
+      // rs2: GridBounds(355,105,512,360)
+      val rasterSources = List(
+        rasterGeoTiffPath("vlm/lc8-utm-1.tif"),
+        rasterGeoTiffPath("vlm/lc8-utm-2.tif")
+      ).map(GeoTiffRasterSource(_))
+
+      val mosaic = MosaicRasterSource.instance(NonEmptyList.fromListUnsafe(rasterSources), rasterSources.head.crs)
+
+      val layout = mosaic
+        .reproject(WebMercator)
+        .tileToLayout(ld, identity, NearestNeighbor, AutoHigherResolution)
+
+      val mosaicReprojected = mosaic
+        .reproject(WebMercator)
+        .resampleToGrid(ld, NearestNeighbor, AutoHigherResolution)
+
+      layout.keys.foreach { key =>
+        val extent = ld.mapTransform.keyToExtent(key)
+        val ltile = layout.read(key).map(_.band(0).toArrayTile)
+        val mtile = mosaicReprojected.read(extent.buffer(mosaicReprojected.cellSize.resolution / 2)).map(_.mapTile(_.band(0).toArrayTile))
+
+        (ltile, mtile) match {
+          case (Some(ltile), Some(mtile)) =>
+            val sample = Raster(ArrayTile.alloc(ltile.cellType, 256, 256), extent)
+            val rtile = sample.merge(mtile)
+            assertTilesEqual(ltile, rtile.tile)
+          case (None, None) =>
+          case (l, r) => throw new Exception(s"$l is not equal to $r")
+        }
+      }
+    }
   }
 }

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

@@ -93,7 +93,31 @@ abstract class MosaicRasterSource extends RasterSource {
   }
 
   def read(bounds: GridBounds[Long], bands: Seq[Int]): Option[Raster[MultibandTile]] = {
-    val rasters = sources map { _.read(bounds, bands) }
+    /** The passed bounds are relative to the [[MosaicRasterSource]] bounds.
+      * However, each [[RasterSource]] has its own [[GridBounds]].
+      * Before passing [[GridBounds]] into each underlying [[RasterSource]]
+      * we need to map them into the each [[RasterSource]] relative grid space.
+      *
+      * This is done by calculating the relative offset using each [[RasterSource]]
+      * underlying [[Extent]].
+
+      * @param gb     global bounds, relative to the [[MosaicRasterSource]]
+      * @param extent extent of the [[RasterSource]]
+      * @return       relative to the extent [[GridBounds]]
+      */
+    def relativeGridBounds(gb: GridBounds[Long], extent: Extent): GridBounds[Long] = {
+      val GridBounds(colMin, rowMin, colMax, rowMax) = gb
+      val (sourceColOffset, sourceRowOffset) = gridExtent.mapToGrid(extent.xmin, extent.ymax)
+
+      GridBounds(
+        colMin - sourceColOffset,
+        rowMin - sourceRowOffset,
+        colMax - sourceColOffset,
+        rowMax - sourceRowOffset
+      )
+    }
+
+    val rasters = sources.map { rs => rs.read(relativeGridBounds(bounds, rs.extent), bands) }
     rasters.reduce
   }