Add union to merge methods

CHANGELOG.md

@@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Added
 - Add RasterSourceRDD.tiledLayerRDD within the geometry intersection [#3474](https://github.com/locationtech/geotrellis/pull/3474)
 - Expose AWS_REQUEST_PAYER environment variable [#3479](https://github.com/locationtech/geotrellis/pull/3479)
+- Add union/outer join merger of tiles [#3482](https://github.com/locationtech/geotrellis/pull/3482)
 
 ### Changed
 - Migration to CE3 and other major dependencies upgrade [#3389](https://github.com/locationtech/geotrellis/pull/3389)

raster/src/main/scala/geotrellis/raster/merge/MultibandTileMergeMethods.scala

@@ -63,4 +63,35 @@ trait MultibandTileMergeMethods extends TileMergeMethods[MultibandTile] {
 
     ArrayMultibandTile(bands)
   }
+
+    /**
+    * Union this [[MultibandTile]] with the other one.  The output tile's
+    * extent will be the minimal extent which encompasses both input extents.
+    * A new MutlibandTile is returned.
+    *
+    * @param   extent        The extent of this MultiBandTile
+    * @param   otherExtent   The extent of the other MultiBandTile
+    * @param   other         The other MultiBandTile
+    * @param   method        The resampling method
+    * @param   unionFunc     The function which decides how values from rasters being combined will be transformed
+    * @return                A new MultiBandTile, the result of the merge
+    */
+  def union(
+    extent: Extent,
+    otherExtent: Extent,
+    other: MultibandTile,
+    method: ResampleMethod,
+    unionFunc: (Option[Double], Option[Double]) => Double
+  ): MultibandTile = {
+    val bands: Seq[Tile] =
+      for {
+        bandIndex <- 0 until self.bandCount
+      } yield {
+        val thisBand = self.band(bandIndex)
+        val thatBand = other.band(bandIndex)
+        thisBand.union(extent, otherExtent, thatBand, method, unionFunc)
+      }
+
+    ArrayMultibandTile(bands)
+  }
 }

raster/src/main/scala/geotrellis/raster/merge/RasterMergeMethods.scala

@@ -50,4 +50,16 @@ abstract class RasterMergeMethods[
     */
   def merge(other: Raster[T]): Raster[T] =
     merge(other, NearestNeighbor)
+
+  /**
+    * Union this [[Raster]] with the other one.  All places in the
+    * present raster that contain NODATA are filled-in with data from
+    * the other raster.  A new Raster is returned.
+    *
+    * @param   other         The other Raster
+    * @param   method        The resampling method
+    * @return                A new Raster, the result of the merge
+    */
+  def union(other: Raster[T], method: ResampleMethod, unionFunc: (Option[Double], Option[Double]) => Double): Raster[T] =
+    Raster(self.tile.union(self.extent, other.extent, other.tile, method, unionFunc), self.extent.combine(other.extent))
 }

raster/src/main/scala/geotrellis/raster/merge/RasterTileFeatureMergeMethods.scala

@@ -32,4 +32,7 @@ abstract class RasterTileFeatureMergeMethods[
 
   def merge(other: TileFeature[Raster[T], D], method: ResampleMethod): TileFeature[Raster[T], D] =
     TileFeature(self.tile.merge(other.tile, method), Semigroup[D].combine(self.data, other.data))
+
+  def union(other: TileFeature[Raster[T], D], method: ResampleMethod, unionFunc: (Option[Double], Option[Double]) => Double): TileFeature[Raster[T], D] =
+    TileFeature(self.tile.union(other.tile, method, unionFunc), Semigroup[D].combine(self.data, other.data))
 }

raster/src/main/scala/geotrellis/raster/merge/SinglebandTileMergeMethods.scala

@@ -99,52 +99,48 @@ trait SinglebandTileMergeMethods extends TileMergeMethods[Tile] {
         val gridBounds = re.gridBoundsFor(sharedExtent)
         val targetCS = CellSize(sharedExtent, gridBounds.width, gridBounds.height)
 
+        val interpolate = Resample(method, other, otherExtent, targetCS)
+
         self.cellType match {
           case BitCellType | ByteCellType | UByteCellType | ShortCellType | UShortCellType | IntCellType  =>
-            val interpolate: (Double, Double) => Int = Resample(method, other, otherExtent, targetCS).resample _
-            // Assume 0 as the transparent value
             cfor(0)(_ < self.rows, _ + 1) { row =>
               cfor(0)(_ < self.cols, _ + 1) { col =>
                 if (self.get(col, row) == 0) {
                   val (x, y) = re.gridToMap(col, row)
-                  val v = interpolate(x, y)
+                  val v = interpolate.resample(x, y)
                   if(isData(v)) {
                     mutableTile.set(col, row, v)
                   }
                 }
               }
             }
           case FloatCellType | DoubleCellType =>
-            val interpolate: (Double, Double) => Double = Resample(method, other, otherExtent, targetCS).resampleDouble _
-            // Assume 0.0 as the transparent value
             cfor(0)(_ < self.rows, _ + 1) { row =>
               cfor(0)(_ < self.cols, _ + 1) { col =>
                 if (self.getDouble(col, row) == 0.0) {
                   val (x, y) = re.gridToMap(col, row)
-                  val v = interpolate(x, y)
+                  val v = interpolate.resampleDouble(x, y)
                   if(isData(v)) {
                     mutableTile.setDouble(col, row, v)
                   }
                 }
               }
             }
           case x if x.isFloatingPoint =>
-            val interpolate: (Double, Double) => Double = Resample(method, other, otherExtent, targetCS).resampleDouble _
             cfor(0)(_ < self.rows, _ + 1) { row =>
               cfor(0)(_ < self.cols, _ + 1) { col =>
                 if (isNoData(self.getDouble(col, row))) {
                   val (x, y) = re.gridToMap(col, row)
-                  mutableTile.setDouble(col, row, interpolate(x, y))
+                  mutableTile.setDouble(col, row, interpolate.resampleDouble(x, y))
                 }
               }
             }
           case _ =>
-            val interpolate: (Double, Double) => Int = Resample(method, other, otherExtent, targetCS).resample _
             cfor(0)(_ < self.rows, _ + 1) { row =>
               cfor(0)(_ < self.cols, _ + 1) { col =>
                 if (isNoData(self.get(col, row))) {
                   val (x, y) = re.gridToMap(col, row)
-                  mutableTile.set(col, row, interpolate(x, y))
+                  mutableTile.set(col, row, interpolate.resample(x, y))
                 }
               }
             }
@@ -154,4 +150,80 @@ trait SinglebandTileMergeMethods extends TileMergeMethods[Tile] {
       case _ =>
         self
     }
+
+  /** Unions this tile with another tile, preserving the non-overlapping
+    * portions of each tile extent.
+    *
+    * This method requires a union function to be provided which decides
+    * how values will be added to the output raster. A `None` represents
+    * one of the two input rasters not having a value at the resampled location.
+    * Two `None` values mean that both tiles are missing values at the
+    * cell in question.
+    */
+  def union(extent: Extent, otherExtent: Extent, other: Tile, method: ResampleMethod, unionFunc: (Option[Double], Option[Double]) => Double): Tile = {
+    val unionInt: (Option[Int], Option[Int]) => Int =
+      (l: Option[Int], r: Option[Int]) => {
+        unionFunc(l.map(_.toDouble), r.map(_.toDouble)).toInt
+      }
+
+    val combinedExtent = otherExtent combine extent
+    val re = RasterExtent(extent, self.cols, self.rows)
+    val gridBounds = re.gridBoundsFor(combinedExtent, false)
+    val targetCS = CellSize(combinedExtent, gridBounds.width, gridBounds.height)
+    val targetRE = RasterExtent(combinedExtent, targetCS)
+    val mutableTile = ArrayTile.empty(self.cellType, targetRE.cols, targetRE.rows)
+
+    val interpolateLeft = Resample(method, self, extent, targetCS)
+    val interpolateRight = Resample(method, other, otherExtent, targetCS)
+
+    self.cellType match {
+      case BitCellType | ByteCellType | UByteCellType | ShortCellType | UShortCellType | IntCellType  =>
+        cfor(0)(_ < targetRE.rows, _ + 1) { row =>
+          cfor(0)(_ < targetRE.cols, _ + 1) { col =>
+            val (x,y) = targetRE.gridToMap(col, row)
+            val l = interpolateLeft.resample(x, y)
+            val r = interpolateRight.resample(x, y)
+            val maybeL = Some(l)
+            val maybeR = Some(r)
+            val v = unionInt(maybeL, maybeR)
+            mutableTile.set(col, row, v)
+          }
+        }
+      case FloatCellType | DoubleCellType =>
+        cfor(0)(_ < targetRE.rows, _ + 1) { row =>
+          cfor(0)(_ < targetRE.cols, _ + 1) { col =>
+            val (x,y) = targetRE.gridToMap(col, row)
+            val l = interpolateLeft.resampleDouble(x, y)
+            val r = interpolateRight.resampleDouble(x, y)
+            val maybeL = Some(l)
+            val maybeR = Some(r)
+            val v = unionFunc(maybeL, maybeR)
+            mutableTile.setDouble(col, row, v)
+          }
+        }
+      case x if x.isFloatingPoint =>
+        cfor(0)(_ < targetRE.rows, _ + 1) { row =>
+          cfor(0)(_ < targetRE.cols, _ + 1) { col =>
+            val (x,y) = targetRE.gridToMap(col, row)
+            val l = interpolateLeft.resampleDouble(x, y)
+            val r = interpolateRight.resampleDouble(x, y)
+            val maybeL = if (isNoData(l)) None else Some(l)
+            val maybeR = if (isNoData(r)) None else Some(r)
+            mutableTile.setDouble(col, row, unionFunc(maybeL, maybeR))
+          }
+        }
+      case _ =>
+        cfor(0)(_ < targetRE.rows, _ + 1) { row =>
+          cfor(0)(_ < targetRE.cols, _ + 1) { col =>
+            val (x,y) = targetRE.gridToMap(col, row)
+            val l = interpolateLeft.resample(x, y)
+            val r = interpolateRight.resample(x, y)
+            val maybeL = if (isNoData(l)) None else Some(l)
+            val maybeR = if (isNoData(r)) None else Some(r)
+            mutableTile.set(col, row, unionInt(maybeL, maybeR))
+          }
+        }
+    }
+    mutableTile
+  }
 }

raster/src/main/scala/geotrellis/raster/merge/TileFeatureMergeMethods.scala

@@ -19,16 +19,27 @@ package geotrellis.raster.merge
 import geotrellis.raster._
 import geotrellis.raster.resample._
 import geotrellis.vector._
+
 import cats.Semigroup
+import cats.syntax.semigroup._ 
 
 
 abstract class TileFeatureMergeMethods[
   T <: CellGrid[Int]: * => TileMergeMethods[T],
   D: Semigroup
 ](val self: TileFeature[T, D]) extends TileMergeMethods[TileFeature[T, D]] {
   def merge(other: TileFeature[T, D], baseCol: Int, baseRow: Int): TileFeature[T, D] =
-    TileFeature(self.tile.merge(other.tile, baseCol, baseRow), Semigroup[D].combine(self.data, other.data))
+    TileFeature(self.tile.merge(other.tile, baseCol, baseRow), self.data combine other.data)
 
   def merge(extent: Extent, otherExtent: Extent, other: TileFeature[T, D], method: ResampleMethod): TileFeature[T, D] =
-    TileFeature(self.tile.merge(extent, otherExtent, other.tile, method), Semigroup[D].combine(self.data, other.data))
+    TileFeature(self.tile.merge(extent, otherExtent, other.tile, method), self.data combine other.data)
+
+  def union(
+    extent: Extent,
+    otherExtent: Extent,
+    other: TileFeature[T, D],
+    method: ResampleMethod,
+    unionFunc: (Option[Double], Option[Double]) => Double
+  ): TileFeature[T, D] =
+    TileFeature(self.tile.union(extent, otherExtent, other.tile, method, unionFunc), self.data combine other.data)
 }

raster/src/main/scala/geotrellis/raster/merge/TileMergeMethods.scala

@@ -78,4 +78,9 @@ trait TileMergeMethods[T] extends MethodExtensions[T] {
   def merge(extent: Extent, otherExtent: Extent, other: T): T =
     merge(extent, otherExtent, other, NearestNeighbor)
 
+
+  def union(extent: Extent, otherExtent: Extent, other: T, method: ResampleMethod, unionFunc: (Option[Double], Option[Double]) => Double): T
+
+  def union(extent: Extent, otherExtent: Extent, other: T, unionFunc: (Option[Double], Option[Double]) => Double): T =
+    union(extent, otherExtent, other, NearestNeighbor, unionFunc)
 }

raster/src/test/scala/geotrellis/raster/merge/TileUnionMethodsSpec.scala

@@ -0,0 +1,80 @@
+/*
+ * Copyright 2016 Azavea
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+package geotrellis.raster.merge
+
+import geotrellis.raster._
+import geotrellis.raster.testkit._
+import geotrellis.raster.resample.NearestNeighbor
+import geotrellis.vector.Extent
+
+import org.scalatest.matchers.should.Matchers
+import org.scalatest.funspec.AnyFunSpec
+
+class TileUnionMethodsSpec extends AnyFunSpec
+    with Matchers
+    with TileBuilders
+    with RasterMatchers {
+  describe("SinglebandTileMergeMethods") {
+
+    it("should union two tiles such that the extent of the output is equal to the minimum extent which covers both") {
+      val cellTypes: Seq[CellType] =
+        Seq(
+          BitCellType,
+          ByteCellType,
+          ByteConstantNoDataCellType,
+          ByteUserDefinedNoDataCellType(1.toByte),
+          UByteCellType,
+          UByteConstantNoDataCellType,
+          UByteUserDefinedNoDataCellType(1.toByte),
+          ShortCellType,
+          ShortConstantNoDataCellType,
+          ShortUserDefinedNoDataCellType(1.toShort),
+          UShortCellType,
+          UShortConstantNoDataCellType,
+          UShortUserDefinedNoDataCellType(1.toShort),
+          IntCellType,
+          IntConstantNoDataCellType,
+          IntUserDefinedNoDataCellType(1),
+          FloatCellType,
+          FloatConstantNoDataCellType,
+          FloatUserDefinedNoDataCellType(1.0f),
+          DoubleCellType,
+          DoubleConstantNoDataCellType,
+          DoubleUserDefinedNoDataCellType(1.0)
+        )
+
+      for(ct <- cellTypes) {
+        val arr = Array.ofDim[Double](100).fill(5.0)
+        arr(50) = 1.0
+        arr(55) = 0.0
+        arr(60) = Double.NaN
+
+        val tile1 =
+          DoubleArrayTile(arr, 10, 10).convert(ct)
+        val e1 = Extent(0, 0, 1, 1)
+        val tile2 =
+          tile1.prototype(ct, tile1.cols, tile1.rows)
+        val e2 = Extent(1, 1, 2, 2)
+        val unioned = tile1.union(e1, e2, tile2, NearestNeighbor, (d1, d2) => d1.getOrElse(4))
+        withClue(s"Failing on cell type $ct: ") {
+          unioned.rows shouldBe (20)
+          unioned.cols shouldBe (20)
+        }
+      }
+    }
+  }
+}