Add union to merge methods

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   unionF      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,
+    unionF: (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, unionF)
+      }
+
+    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, unionF: (Option[Double], Option[Double]) => Double): Raster[T] =
+    Raster(self.tile.union(self.extent, other.extent, other.tile, method, unionF), 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, unionF: (Option[Double], Option[Double]) => Double): TileFeature[Raster[T], D] =
+    TileFeature(self.tile.union(other.tile, method, unionF), Semigroup[D].combine(self.data, other.data))
 }

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

@@ -154,4 +154,69 @@ trait SinglebandTileMergeMethods extends TileMergeMethods[Tile] {
       case _ =>
         self
     }
+
+  def union(extent: Extent, otherExtent: Extent, other: Tile, method: ResampleMethod, unionF: (Option[Double], Option[Double]) => Double): Tile = {
+    val unionInt = (l: Option[Double], r: Option[Double]) => unionF(l, r).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)
+
+    self.cellType match {
+      case BitCellType | ByteCellType | UByteCellType | ShortCellType | UShortCellType | IntCellType  =>
+        val interpolateLeft: (Double, Double) => Int = Resample(method, self, extent, targetCS).resample _
+        val interpolateRight: (Double, Double) => Int = Resample(method, other, otherExtent, targetCS).resample _
+        // Assume 0 as the transparent value
+        cfor(0)(_ < targetRE.rows, _ + 1) { row =>
+          cfor(0)(_ < targetRE.cols, _ + 1) { col =>
+            val (x,y) = targetRE.gridToMap(col, row)
+            val (l,r) = (interpolateLeft(x, y), interpolateRight(x, y))
+            mutableTile.set(col, row, unionInt(Some(l.toDouble), Some(r.toDouble)))
+          }
+        }
+      case FloatCellType | DoubleCellType =>
+        val interpolateLeft: (Double, Double) => Double = Resample(method, self, extent, targetCS).resampleDouble _
+        val interpolateRight: (Double, Double) => Double = Resample(method, other, otherExtent, targetCS).resampleDouble _
+
+        // Assume 0.0 as the transparent value
+        cfor(0)(_ < targetRE.rows, _ + 1) { row =>
+          cfor(0)(_ < targetRE.cols, _ + 1) { col =>
+            val (x,y) = targetRE.gridToMap(col, row)
+            val (l,r) = (interpolateLeft(x, y), interpolateRight(x, y))
+            mutableTile.setDouble(col, row, unionF(Some(l), Some(r)))
+          }
+        }
+      case x if x.isFloatingPoint =>
+        val interpolateLeft: (Double, Double) => Double = Resample(method, self, extent, targetCS).resampleDouble _
+        val interpolateRight: (Double, Double) => Double = Resample(method, other, otherExtent, targetCS).resampleDouble _
+        cfor(0)(_ < targetRE.rows, _ + 1) { row =>
+          cfor(0)(_ < targetRE.cols, _ + 1) { col =>
+            val (x,y) = targetRE.gridToMap(col, row)
+            val l = interpolateLeft(x, y)
+            val r = interpolateRight(x, y)
+            val maybeL = if (isNoData(l)) None else Some(l)
+            val maybeR = if (isNoData(r)) None else Some(r)
+            mutableTile.setDouble(col, row, unionF(maybeL, maybeR))
+          }
+        }
+      case _ =>
+        val interpolateLeft: (Double, Double) => Int = Resample(method, self, extent, targetCS).resample _
+        val interpolateRight: (Double, Double) => Int = Resample(method, other, otherExtent, targetCS).resample _
+        cfor(0)(_ < targetRE.rows, _ + 1) { row =>
+          cfor(0)(_ < targetRE.cols, _ + 1) { col =>
+            val (x,y) = targetRE.gridToMap(col, row)
+            val l = interpolateLeft(x, y)
+            val r = interpolateRight(x, y)
+            val maybeL = if (isNoData(l)) None else Some(l.toDouble)
+            val maybeR = if (isNoData(r)) None else Some(r.toDouble)
+            mutableTile.set(col, row, unionInt(maybeL, maybeR))
+            //if (l!=r) println(s"x => ${x}, y => ${y}, col => ${col}, row => ${row} | l,r => ${l}, ${r}")
+          }
+        }
+    }
+    mutableTile
+  }
 }

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

@@ -31,4 +31,13 @@ abstract class TileFeatureMergeMethods[
 
   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))
+
+  def union(
+    extent: Extent,
+    otherExtent: Extent,
+    other: TileFeature[T, D],
+    method: ResampleMethod,
+    unionF: (Option[Double], Option[Double]) => Double
+  ): TileFeature[T, D] =
+    TileFeature(self.tile.union(extent, otherExtent, other.tile, method, unionF), Semigroup[D].combine(self.data, 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, unionF: (Option[Double], Option[Double]) => Double): T
+
+  def union(extent: Extent, otherExtent: Extent, other: T, unionF: (Option[Double], Option[Double]) => Double): T =
+    union(extent, otherExtent, other, NearestNeighbor, unionF)
 }

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)
+        }
+      }
+    }
+  }
+}